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ABSTRACT 


NAVIGATION AND CONTROL STUDIES ON CRUISE 
MISSILES 


EKÜTEKİN, Vedat 
Ph. D., Department of Mechanical Engineering 
Supervisor: Prof. Dr. M. Kemal ÖZGÖREN 


January 2007, 298 pages 


A cruise missile is a guided missile that uses a lifting wing and a jet 
propulsion system to allow sustained flight. Cruise missiles are, in essence, 
unmanned aircraft and they are generally designed to carry a large conventional or 
nuclear warhead many hundreds of miles with excellent accuracy. In this study, 
navigation and control studies on cruise missiles are performed. Due to the variety 
and complexity of the subsystems of the cruise missiles, the main concern is limited 
with the navigation system. Navigation system determines the position, velocity, 
attitude and time solutions of the missile. Therefore, it can be concluded that an 
accurate self-contained navigation system directly influences the success of the 
missile. In the study, modern radar data association algorithms are implemented as 
new Terrain Aided Navigation (TAN) algorithms which can be used with low-cost 
Inertial Measurement Units (IMU’s). In order to perform the study, first a thorough 


survey of the literature on mid-course navigation of cruise missiles is performed. 


iv 


Then, study on modern radar data association algorithms and their implementations 
to TAN are done with simple simulations. At the case study part, a six degree of 
freedom (6 DOF) flight simulation tool is developed which includes the 
aerodynamic and dynamic model of the cruise missile model including error model 
of the navigation system. Finally, the performances of the designed navigation 
systems with the implemented TAN algorithms are examined in detail with the help 


of the simulations performed. 


Keywords: Cruise Missile, Terrain Aided Navigation (TAN), Probabilistic Data 
Association Filter (PDAF), Track Splitting Filter (TSF), Multiple Hypothesis 
Tracking (MHT). 


ÖZ 


SEYİR FÜZELERİ ÜZERİNE SEYRÜSEFER VE DENETİM 
ÇALIŞMALARI 


EKÜTEKİN, Vedat 
Doktora, Makina Mühendisliği Bölümü 
Tez Yöneticisi: Prof. Dr. M. Kemal ÖZGÖREN 


Ocak 2007, 298 sayfa 


Seyir füzesi, kaldırma kanatları ve jet itki sistemi ile kararlı uçuş sağlayan 
güdümlü bir füzedir. Seyir füzeleri genellikle, büyük konvansiyonel ya da nükleer 
savaş başlıklarını uzak mesafelere çok hassas olarak taşıyan insansız hava 
taşıtlarıdır. Bu çalışmada, seyir füzeleri üzerine seyrüsefer ve denetim çalışmaları 
gerçekleştirilmiştir. Seyir füzelerinin alt sistemlerindeki çeşitlilik ve karmaşıklık 
nedeniyle, çalışmanın ana konusu seyrüsefer sistemiyle kısıtlanmıştır. Seyrüsefer 
sistemi füzenin konum, hız, yönelim ve zaman çözümlerini belirler. Bu nedenle, 
hassas, kendi kendine yeterli bir seyrüsefer sisteminin füzenin başarısını doğrudan 
etkileyeceği sonucuna varılabilir. Çalışmada, modern radar veri ilişkilendirme 
algoritmaları, düşük maliyetli ataletsel seyrüsefer sistemleri ile kullanılabilecek yeni 
Arazi Destekli Seyrüsefer (ADS) algoritma uygulamaları için kullanılmıştır. 
Çalışmayı gerçekleştirmek için, ilk aşamada seyir füzelerinin seyrüsefer 


yöntemlerine ait ayrıntılı kaynak araştırması yapılmıştır. Daha sonra, modem radar 


vi 


veri ilişkilendirme algoritmaları üzerine çalışılmış ve bunların ADS uygulamaları 
basit benzetimlerle gerçekleştirilmiştir. Örnek olay incelemesi kısmında, seyir 
füzesinin aerodinamik ve dinamik modeli ile seyrüsefer sisteminin hata 
modellemesini de içeren altı serbestlik dereceli bir uçuş benzetim aracı 
geliştirilmiştir. Son olarak, yeni uygulanan ADS algoritmaları kullanılarak 
tasarlanan seyrüsefer sistemlerinin — başarımları, gerçekleştirilen uçuş 


benzetimlerinin yardımıyla ayrıntılı olarak incelenmiştir. 


Anahtar Kelimeler: Seyir Füzesi, Arazi Destekli Seyrüsefer (ADS), Olasılıklı Veri 
İlişkilendirme Filtresi (OVIF), İz Ayırma Filtresi (İAF), Çoklu Varsayımlı Takip 


(ÇVT) 
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CHAPTER 1 


INTRODUCTION 


1.1. Scopeof the Study 


A cruise missile is a guided missile that uses a lifting wing and a jet 
propulsion system to allow sustained flight. Cruise missiles are, in essence, 
unmanned aircraft and they are generally designed to carry a large conventional or 


nuclear warhead many hundreds of miles with excellent accuracy [1]. 


A cruise missile usually flies at subsonic speed and it would reguire several 
hours of continuously guided flight to cover its mission distance. Hence, guidance 
errors that accumulate with time would be almost 100 times larger for a cruise 
missile than for a ballistic missile which is guided for the first five of the twenty 
minutes. Therefore its accurate arrival on the target could be achieved only with 
continuous guidance that is updated and corrected from time to time by new 
location information. In order to obtain the necessary location information, a long- 
range cruise missile employs a device that can correlate information obtained by an 
onboard sensor about the terrain it is flying over with some kind of map stored in 


the memory of an onboard computer [2]. 


Navigation system of a cruise missile determines the position, velocity, 
attitude and time solutions of the missile. Therefore, it can be concluded that an 
accurate self-contained navigation system directly influences the success of the 
missile. In this study, navigation and control studies on cruise missiles will be 


performed. 


Terrain Aided (Referenced) Navigation (TAN) is an important part of 
“Integrated Navigation Systems” in military and civil avionics. TAN provides 
position fixes, which can be used to aid a central navigation system. Especially, if 
other sources for position aids, like the Global Positioning System (GPS), are not 
available, TAN can provide reliable position information in low level flights over 


significant terrain [3]. 


The scope of the study is to implement some modern radar data association 
algorithms as new Terrain Aided Navigation (TAN) algorithms which can be used 


with low-cost Inertial Measurement Units (IMU’s). 


In this chapter, theory about the study will be given. First, information about 
cruise missiles and cruise missile navigation performance will be given. Then, 
literature survey on TAN techniques will be discussed in detail. Finally, information 
about radar tracking techniques and possible implementations of radar data 
association algorithms to TAN will be given. At the last section of the chapter, 


outline of the thesis study will be summarized. 


12. Cruise Missiles 


1.2.1. Background 


A cruise missile is a guided missile that uses a lifting wing and a jet 
propulsion system to allow sustained flight. Cruise missiles are, in essence, 
unmanned aircraft. They are generally designed to carry a large conventional or 
nuclear warhead many hundreds of miles with excellent accuracy. In 2001, modern 
cruise missiles normally travel at sub-sonic speeds, are self-navigating, and fly low 
in order to avoid radar detection [1]. The term cruise missile covers several vehicles 
and their capabilities, from the Chinese Silkworm (HY-2), which has a range of less 
than 105 km, to the U.S. Advanced Cruise Missile (ACM), which can fly to ranges 
of up to 3,000 km. These vehicles vary greatly in their speed and ability to penetrate 
defenses. All, however, meet the definition of a cruise missile: “an unmanned self- 
propelled guided vehicle that sustains flight through aerodynamic lift for most of its 
flight path and whose primary mission is to place an ordnance or special payload on 
a target”. This definition can include unmanned air vehicles (UAV’s) and 


unmanned control-guided helicopters or aircraft [4]. 


Cruise missiles were first developed by Nazi Germany during World War II. 
The V-1 (introduced in 1944) was the first weapon to use the classic cruise missile 
layout of a bomb-like fuselage with short wings and a dorsally mounted engine, 
along with a simple inertial guidance system. The V-1 was propelled by a crude 
pulse-jet engine, the sound of which gave the V-1 its nickname of “buzz bomb”. 
Japanese kamikaze aircraft could be viewed as manned cruise missiles. During the 
Cold War, both the United States and the Soviet Union experimented further with 
the concept, deploying early cruise missiles from submarines and aircraft. The 
Soviet Union was especially fond of large cruise missiles. The United States had a 
program to develop a nuclear-powered cruise missile, Project Pluto. Although the 
concept was proven sound, none were ever test-launched. While ballistic missiles 


were the weapons of choice for land targets, heavy nuclear and conventional tipped 


cruise missiles were seen by the USSR as a primary weapon to destroy US carrier 
battle groups. Large submarines (e.g. Echo and Oscar class) were developed to 
carry these weapons and shadow US battle groups at sea, and large bombers (e.g. 


Backfire, Bear, and Blackjack models) were eguipped with the weapons [1]. 


As of 2001, the Tomahawk missile (BGM-109) model has become a 
significant part of the US naval arsenal. It gives ships and submarines an extremely 
accurate, long-range, conventional land attack weapon. Each costs about $1,000,000 
USD. The United States Air Force deploys an air launched cruise missile, the 
AGM-86. It can be launched from bombers like the B-52 Stratofortress. Both the 


Tomahawk and the AGM-86 were used extensively during Operation Desert Storm 


[1]. 


In Figure 1, Figure 2, Table 1 and Table 2, well-known US cruise missiles 


with specifications were presented as examples to typical cruise missiles. 








Figure 2. AGM-86 Air-Launched Cruise Missile (ALCM) (6) 


1.2.2. Cruise Missile Technology 


Cruise missile technology has advanced substantially since the German V-1 
of World War II. Modern cruise missiles fly at altitudes one-tenth those of the V-1, 
have Radar Cross-Sections (RCS) one hundred times smaller (which reduces 


detectability), and accuracies two hundred times better [7]. 


The technology of the cruise missile has four main component elements: 


1. Airframe; 


2. Propulsion system; 


3. Guidance systems; 


4. Warhead. [7] 


As an example, the main components of Tomahawk cruise missile are 


shown in Figure 3 [5]. 


Table 1. BGM-109 Tomahawk Cruise Missile Specifications [5] 








Primary Function: 


Long-range subsonic cruise missile for attacking land targets. 





Contractor: 


Hughes Missile Systems Co., Tucson, Ariz. 





Power Plant: 


Williams International F107-WR-402 cruise turbo-fan engine; 
solid-fuel booster 





18 feet 3 inches (5.56 meters); with booster: 20 feet 6 inches 














en (6.25 meters) 
Weight: 2,650 pounds (1192.5 kg); 3,200 pounds (1440 kg) with 
booster 

Diameter: 20.4 inches (51.81 cm) 

Wing Span: 8 feet 9 inches (2.67 meters) 

R ; Land attack, conventional warhead: 600 nautical miles (690 
EN statute miles, 1104 km) 

Speed: Subsonic — about 550 mph (880 km/h) 


Guidance System: 


Inertial and TERCOM 





Conventional: 1,000 pounds Bull pup, or conventional sub- 
munitions dispenser with combined effect bomblets, or WDU- 








TS 36 warhead with PBXN-107 explosive & FMU-148 fuze, or 
200 kt. W-80 nuclear device 

Date Deployed: 1983 
$500,000 - current production Unit Cost 

Costs $1,400,000 - average unit cost (TY$) 


$11,210,000,000 - total program cost (TY$) 











Total Program 





4 170 missiles 














Table 2. AGM-86 Air-Launched Cruise Missile (ALCM) Specifications [6] 








Primary Function: 


Air-to-surface strategic missile 








Contractor: Boeing Aerospace Co. 
Guidance Litton Guidance and Control 
Contractors: 





Power Plant: 


Williams Research Corp. F-107-WR-10 turbofan engine 




















Thrust: 600 pounds (270 kilograms) 

Length: 20 feet, 9 inches (6.29 meters) 

Weight: 3,150 pounds (1,417.5 kilograms) 

Diameter: 24.5 inches (62.23 centimeter) 

Wingspan: 12 feet (3.64 meters) 

Range: AGM-86B: 1,500-plus miles (1,305 nautical miles) 
Speed: About 550 mph (Mach 0.73) 





Guidance System: 


Litton inertial navigation element with terrain contour- 
matching updates 























Warheads: Nuclear capable 
A terrain contour-matching guidance system that allows the 
Sensors: missile to fly complicated routes to a target through use of 
maps of the planned flight route stored in on-board computers 
Unit Cost: $1 million 
Date Deployed: December 1982 
Inventory: Active force, 1,628; ANG, 0; Reserve, 0 
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Figure 3. Main Components of Tomahawk Cruise Missile 


1.2.2.1. Airframe 


The airframe is essentially that of a small (unmanned) aeroplane or a design 
based on a Remotely Piloted Vehicle (RPV). Early two-wing, three-surface tail 
aircraft designs were followed by four-wing, four-tail cruciform configurations. As 
an example, the body of a Tomahawk cruise missile which has two wings and a 


four-fan tail [7] is given in Figure 4 [5]. 





Figure 4. Tomahawk Cruise Missile Configuration 


1.2.2.2. Propulsion System 


The propulsion system needs to maintain sufficient momentum to counter 
the force of gravity. Most cruise missiles are propelled by a small, highly 
specialized, air-breathing engine which thus needs to draw oxygen from the 
atmosphere into the engine for the bulk of the flight. Air-breathing engines are one 
of four types: pulsejet, ramjet, turbojet or turbofan, this last being a more efficient 
form of turbojet developed in the 1970’s. Short to medium range systems tend to 
employ turbojets, which though less efficient, are usually less expensive than 
turbofans. Most long-range missiles, e.g. the Tomahawk use highly efficient 
turbofans propelling them at high subsonic speeds. The few long-range cruise 
missiles propelled by ramjets include the French ASMP and ASURA which are 
capable of Mach 2 and Mach 3 speeds respectively [7]. 


Many missiles are launched by rocket boosters and some missiles, especially 
short-range Anti Ship Cruise Missiles (ASCM’s) like the Exocet, are powered 
throughout their flight by rocket motors. Older rocket-propelled models, such as the 
Styx and Silkworm ASCM’s use liquid fuelled rocket engines, while newer ones, 


such as the Exocet, use solid fuel motors [7]. 


1.2.2.3. Guidance Systems 


Cruise missiles have at least two guidance systems: an in-flight guidance 
system to maintain its flight path and altitude, and a terminal guidance system for 
the final approach to the target. Depending upon the particular characteristics of the 
guidance system, the missile may be programmed: as autonomous (i.e. launch and 
leave); or for remote piloting by command (i.e. flow by a human operator over a 
remote communications link); or as semi-autonomous (a combination of the two, 


with remote manual input in the terminal stage) [7]. 


1.2.2.3.1. In-flight Guidance 


In-flight guidance relies on Inertial Navigation Systems (INS) using 
gyroscopes to ascertain the missile's position. Shorter-range cruise missiles may use 
only inertial and terminal guidance. Longer-range missiles reguire supplemental 
information to make up for inherent inertial guidance inaccuracies (or drift). One 
sophisticated supplemental system in current use is Terrain Contour Matching, 
known in the United States as TERCOM, a position fixing technigue. A digital 
terrain map of the missile’s planned route has first to be made, and it is then stored 
in the weapon's guidance system. Updates received from a radar altimeter 
determine the missile’s altitude and this information is then compared with terrain 
heights in the pre-stored digital map. Once the updates are received, the missile can 


correct its flight back to the planned route [7]. 


Due to the high cost and complexities of obtaining the satellite data needed 
to create the digitized maps for TERCOM, the US is the only nation that currently 
incorporates this technology widely in its cruise missiles. However, the Soviet SS- 
N-21 Sampson, a long-range Submarine-Launched Cruise Missile (SLCM) dubbed 
*Tomahawkski” due to its similarity to the US system, is believed to be able to 
incorporate a TERCOM-like guidance system. The French ASMP and Apache 
missiles also use terrain matching and several other Western nations, including the 
UK and Sweden, are believed to have the capability to incorporate TERCOM into 


cruise missiles [7]. 


Global Positioning System (GPS) data is another supplementary guidance 
system, which has yet to be fully exploited. It uses a constellation of 24 
continuously transmitting navigation satellites provided by the US Department of 
Defense. Military users can receive positional data accurate down to 5 meters or 
less. However, the same signal ‘degraded’ for civilian users (and potential 
adversaries) is less precise, but accurate enough for most purposes. A similar 


system, GLONASS, is being deployed by Russia [7]. 
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GPS is now used in US cruise missiles, such as the Tomahawks used to 
attack Bosnian Serb targets in September 1995. France intends to use GPS for its 
Apache series of cruise missiles, and other countries are also expected to do so. 
Thus, with the ever-widening availability of this technology, all proliferators can 
significantly enhance the accuracy of their cruise missiles. For them GPS will be 
preferable to TERCOM because it does not reguire such elaborate and cost- 


intensive pre-programming of data [7]. 


However both systems have their limitations. TERCOM navigates by 
identifying distinctive terrain features. On the other hand, the GPS system utilizes 
long-range satellite systems whose transmissions can be jammed using shorter 
range, more powerful signals, and the civilian signals can also be switched off if 


necessary, as was the case during Operation Desert Storm in 1991 [7]. 


1.2.2.3.2. Terminal Guidance 


Terminal guidance systems help the missile to home in on the target in the 
final stages of flight. These systems may make use of active or semi-active radar, 


infrared, television, or “home-on-jam” (i.e. on a jamming signal) techniques [7]. 


The Tomahawk uses an additional set of precise terminal navigation updates 
known as the Digital Scene Matching Area Correlator (DSMAC), a two- 
dimensional, map-matching concept that employs an onboard sensor to obtain a 
seguence of images of the ground directly below the missile. The images are 
compared to reference data stored in the missile’s navigational computer, and route 


changes are made accordingly, prior to final target acguisition [7]. 


In Figure 5, TERCOM and DSMAC guidance principles are shown [5]. 
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Figure 5. TERCOM and DSMAC Guidance 


1.2.2.4. Warhead 


In some respects the most significant component of a cruise missile is its 
warhead. Unlike a ballistic missile, which places enormous stresses on its warhead 
as it accelerates and as it re-enters the Earth’s atmosphere, a cruise missile flies 
much like an aircraft. Its warhead can therefore be based upon munitions originally 
designed for manned aircraft, making the development of Chemical and Biological 


(CW and BW) payloads for cruise missiles a comparatively simple matter [7]. 


Currently, most cruise missiles are armed with conventional, high-explosive 
warheads. However several countries are known to have available blast 


fragmentation warheads for use with their cruise missiles (e.g. ASCM’s such as the 
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Chinese Silkworm, the Iragi FAW series, and the Israeli Gabriel) or sub-munitions 
(e.g. the French Exocet ASCM and Apache TLACM and the German Kormoran 
ASCM) [7]. 


Some ex-Soviet systems are dual-capable and can be fitted with either a 
conventional or nuclear warhead. Only the US, Russia and France are known to 
deploy nuclear-armed cruise missiles at present. The Chinese were also reported 
some time ago to be nearing completion of a nuclear warhead for their Silkworm 
ASCM. No nation is currently known to possess a CW or BW cruise missile 
warhead. However there have been media reports suggesting that Syria, Iran, and 


China are attempting to develop these [7]. 


Different variants of the Tomahawk may be nuclear or conventionally 
armed. Due partly to arms control constraints and partly to improvements in 
conventional payloads, the US is concentrating on the development and deployment 
of conventionally armed missiles. The Tomahawk TLAM-D for example carries a 
sub-munitions dispenser that allows it to deliver bomblets on three different targets, 


before diving into a fourth [7]. 


1.2.33. Low-Cost Cruise Missiles 


Advances in new commercial technologies make the development of low 
cost guided weapons possible. US authorities developed the Low Cost Cruise 
Missile Defense (LCCMD) program in order to defeat a threat consisting of 
unsophisticated air vehicles attempting to overwhelm their defensives by attacking 


in large numbers or by attacking over wide geographic areas [8]. 
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It is claimed that 82 countries (including third world countries like Egypt, 
Chile, and Singapore) possess cruise missiles where 75 systems are in service and 


42 are in development (9). 


Advanced low cost interceptor seekers, using commercial hardware, and 
matching seeker performance as Noise Radar Seeker, the Micro Electromechanical 
Machine System (MEMS) Electronically Steerable Antenna (ESA) Seeker, Laser 
Seeker, Infrared (IR) Seeker, Optical ESA, Ultra High Freguency (UHF) Seeker and 
advanced navigation algorithms can be examples for low cost systems which can be 


used in cruise missiles in the future. 


As an example, in reference [10], an amateur researcher claims to build a 
cruise missile in his own garage with a budget of just 5,000 US dollars. He also 


subscribes all the work he does in his site. 


1.3. Terrain Aided Navigation (TAN) 


Terrain Aided Navigation (TAN) is a technique to estimate the position of a 
moving vehicle by comparing the measured terrain profile under the vehicle to a 
stored elevation map. TAN has been operational for unmanned vehicles for some 
time. Although this operational system has proven to be reliable and cost effective, 
it is desirable to develop enhancements which can either reduce the pre-planning 
effort or increase the operational envelope, i.e., reliable operation in terrain with 
less reliable or with stored elevation data with larger errors. It is anticipated that 
terrain aided navigation will be in use for many years to come due to the long term 
stability of the terrain profile of earth, the relative ease of mapping and maintaining 
maps of large operational areas, the ease and reliability with which on-board 
measurements can be made and the relatively low computational burden of 


computing navigation updates in an embedded vehicle processor [11]. 
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TAN is an important part of “Integrated Navigation Systems” in military and 
civil avionics. TAN provides position fixes, which can be used to aid a central 
navigation system. Especially, if other sources for position aids, like the Global 
Positioning System (GPS), are not available, TAN can provide reliable position 
information in low level flights over significant terrain. The outage of the GPS in 
hostile jammed environments or due to shadowing effects caused by low level 
flights in valleys is always possible and has to be expected. Therefore, TAN, which 
is independent from external information sources, is predestinated for additional 


position aiding [3]. 


1.3.1. TAN Technigues 


A number of TAN technigues have been developed and tested. These fall 


into two general algorithmic categories [11]: 


1. Batch Algorithms, 


2. Recursive Algorithms. 


In addition, there are two general map storage technigues: small, high 
fidelity maps which are used at specific points along the intended route of the 
vehicle; and a single, large, low fidelity map which encompasses the entire 
operating area of the vehicle. These technigues are shown in Figure 6 and are 


associated with the two most widely understood TAN implementations [11]: 


1. Terrain Contour Matching (TERCOM) 


2. Sandia Inertial Terrain Aided Navigation (SITAN) [12] 


15 


The most widely known form of TAN is TERCOM. With TERCOM a strip 
of terrain elevation measurements are collected while the vehicle flies along the 
intended route and the measurements are post processed by a batch algorithm to 
provide a correlation with a high fidelity map. In the operational missile systems 
employing TERCOM the stored map preparation and validation process includes 
extensive analysis to evaluate the probability of obtaining a strong and 
unambiguous correlation with candidate maps. The map size in the cross-track 
direction is determined by the “worst case” navigation uncertainty and in the down- 
track direction by the larger of “worst case” navigation uncertainty or the map 
length necessary to provide an unambiguous update opportunity. A seguence of 
maps are then developed to provide navigation update opportunities from the launch 
point to the target. The operational TERCOM applications use a mean absolute 
difference (MAD) algorithm which is only a modest computational reguirement in 
an embedded flight processor. In addition, the map storage reguirements are 
minimized by carefully selecting the minimum number and size of maps reguired 


for each mission [11]. 


In the late 1970’s TAN in the form of SITAN was proposed. SITAN uses an 
extended Kalman Filter (EKF) and a local terrain linearization technigue to 
implement a recursive algorithm. This algorithm operates on individual terrain 
elevation measurements as they become available and for the entire duration of the 
mission. This reguires a map for the entire mission. For missile applications the 
map could be for the length of the mission with the width determined by navigation 
uncertainly and terrain unigueness or suitability. However, for manned aircraft 
applications map data for the entire operating area must be stored because the pilot 
may deviate from the preplanned route at any time. SITAN has been developed and 
evaluated for the manned aircraft application using Digital Terrain Elevation Data 
(DTED) which is a low fidelity Defense Mapping Agency (DMA) product readily 


available in most operational areas [11]. 
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Figure 6. Terrain Aided Navigation (TAN) Technigues (11) 


Other TAN technigues including TERPROM and SPARTAN have been 
developed and evaluated since then. In the late 1950’s and throughout the 1960’s 
when TAN concepts were originally developed and in the 1970’s when TAN 
concepts were applied to missile applications, digital computer capabilities were 
limited. Within the past few years the computational, data storage and memory 
access capabilities of embedded vehicle computers have improved dramatically. 
Thus, the previously assumed computational constraints do not apply as technigues 
are developed to enhance the performance and to expand the operational envelope 


of TAN technigues [11]. 


The TERCOM and SITAN approaches both have attributes that are of 
interest. Although enhancements can be envisioned in a number of areas the 
approach here is to investigate algorithm technigues which would make more 
complete use of the information content of the stored elevation data, the a priori 


knowledge of the errors in the stored elevation data and the elevation 
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measurements. Batch processing algorithm (i.e. TERCOM) was generally selected 
for cruise missiles because: no linearization of the terrain profile is necessary; it 
does not reguire an acguisition process; and the algorithm technigues are applicable 
to both small discrete maps and large maps which can support continuous 


navigation updating [11]. 


1.3.2. TAN System Considerations 


In terrain aided navigation, position estimates are referenced to the terrain 
data and are insensitive to position bias errors in the terrain data. Because of this 
characteristics, terrain aided navigation systems are especially useful in applications 
that reguire accurate navigation relative to targets, obstacles, structures, and other 
features whose locations are derived from the same source as the stored elevation 


data [13]. 


Terrain aided navigation (TAN) consists of sensing a terrain elevation 
profile beneath an air vehicle and correlating the profile with stored digital terrain 
elevation data (DTED) to produce an estimate of vehicle position. An INS, usually 
with barometric altimeter aiding, provides the approximate trajectory. TAN systems 
provide three dimensional position updates to the navigation system by estimating 
INS trajectory errors. Radar or laser altimeter measures ground clearance and the 
DTED gives terrain elevation above mean sea level (MSL). Implementation 
requires an INS, an altimeter, DTED, and a flight computer for executing the TAN 


algorithm. In Figure 7, an illustration is given for TAN measurement process [13]. 
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Figure 7. TAN Measurements 


1.3.3. Digital Terrain Elevation Data (DTED) 


As it is mentioned in the previous sections, the critical part of TAN is the 
elevation model used in the system. Generally, DTED is used for military purposes. 
The U. S. Department of Defense, through the National Geospatial Intelligence 
Agency, produces several kinds of digital cartographic data. One is digital elevation 
data, in a series called DTED. The data is available as 1-by-1 degree guadrangles at 
horizontal resolutions ranging from about 1 kilometer to 1 meter. The lowest 


resolution data is available to the public [14]. 
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DTED Level 0 files have 121-by-121 points. DTED Level 1 files have 1201- 
by-1201. The edges of adjacent tiles have redundant records. DTED files are binary. 
No line ending conversion or byte-swapping is reguired when downloading a 
DTED file [14]. The data available to the public is called Level 0 and has a 30 arc 
second spacing. Other higher resolution data called Level 1 and Level 2 is not 
available to public. Performance specifications of DTED files are defined in a US 


military standard [15] and detailed information can be obtained from there. 


For cruise missile mid-course navigation phase, because of its broad-area 
coverage, Level 1 DTED is used by most TAN systems. With very accurate and 
expensive-to-produce DTED, TAN system horizontal position accuracies rivaling 
those of GPS can be achieved. In TAN systems using Level 1 DTED over broad 
areas, accuracies in the range of 50-200 m CEP are typical for low-flying air 


vehicles like cruise missiles [13]. 


In Turkey, various levels of DTED are prepared by HGK (Harita Genel 
Komutanlığı - Turkish General Mapping Commandership) for all regions of Turkey 
from topographic maps and they are served to national institutions with protocols. 


The properties of DTED prepared for Turkey are given in Table 3 [16]. 


1.4. Literature Survey on TAN 


TAN systems are generally used for military purposes. As a result of this, 
access to literature about TAN became very difficult. Especially for TERCOM, 
original famous report of Baker and Clem (1977), named “Terrain Contour 
Matching (TERCOM) Primer” could not be obtained. However, all the papers about 
TAN found in IEEE and AIAA are investigated and classified for the study. 


Moreover, US patents about TAN are also investigated. 
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Table 3. Properties of DTED Prepared by HGK 



































DTED Type DTED Level 2 DTED Level 1 
Map Scale 1/25,000 1/250,000 
Map Datum WGS84, ED50 WGS84, ED50 
Map Coverage 17x1” 37x3” 

Unit Map Coverage TRS? 19x19 

Unit File Size 0.5 MB 3 MB 
Resolution and +26 m horizontal +130 m horizontal 
Accuracy +20 m altitude +30 m altitude 
DTED Preparation YUKPAF25 YUKPAF250 
Source and Method Interpolation Interpolation 
Confidentiality Classified Unclassified 




















TAN papers can be classified according to their subjects as follows: 


1. Cruise Missile System Performance, 


2. Terrain Models and Path Optimization, 


3. TAN Applications, 


4. TAN Algorithms. 


As it can be seen from the survey results, TAN can be found in various 
subjects related with cruise missiles and navigation applications. The papers are 


investigated considering TAN point of view for the study. 


21 


1.4.1. Cruise Missile System Performance 


In several papers, general system performances of TAN applications of 
cruise missiles are investigated. Henley [17] provides an overview of the 
SPARTAN technigue and other technigues for improving navigation performance 
over very flat terrains. Navigation system performances of various TAN systems are 
investigated and terrain data reguirements are defined. Details of the TAN 
algorithms are not given in the paper. In Table 4, navigation performances of 


various INS aided systems are compared. 


Table 4. Navigation System Performances of Various INS Aiding Systems [17] 
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Nielson [18] investigates the Conventional Air Launched Cruise Missile 
(CALCM) performance. Advantages of integrating GPS navigation into the missile 
in place of TERCOM are stated in this paper. Results and benefits of the GPS 
integrated cruise missile are given considering the applications in the Gulf War. 
However, it is known that jamming is a very important problem for GPS integrated 
systems. In Iraq War, several GPS aided cruise missiles have been jammed by Iraq 
military forces. Therefore, besides ease of using GPS, reliability problems should 


also be taken into consideration. 
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Hicks [19] provides a functional description of the navigation and guidance 
system in the Advanced Cruise Missile (ACM) and discusses some of the areas of 
improvements over the ALCM. From the paper, it can be seen that ACM has a very 
complex navigation and flight control system. The paper is helpful for 


understanding the navigation and guidance of a cruise missile. 


Bennett [20] investigates the use of digital terrain map data for airborne 
operations. The fundamental uses of digital map data for TAN and simulators are 
given in the paper. Moreover, information about mission planning and simulation is 
also given. From his work, Bennett [20] concludes that GPS and TAN are 
complementary navigation sensors, and when properly integrated with INS, they 
provide the essential correlation of aircraft position with respect to the actual 


ground contours. 


1.4.2. Terrain Models and Path Optimization 


Terrain models used in TAN systems are very critical. They should be 
modeled as accurate as possible in order to obtain better navigation solutions. Chen 
and Yu [21] improve the models used for TAN. Actually, they improve the terrain 
model by considering horizontal position noises of the terrain model as colored 
noises. The TAN algorithm used in the paper is SITAN which will be discussed in 


the following sections in detail. 


A very similar paper is presented by Wang and Chen [22]. In their paper, 
possible error sources related to the elevation model is added to the SITAN 
eguations in order to improve navigation solutions. Yu, et al [23] also propose 
various terrain linearization technigues reguired for SITAN implementation and 


present the improvements in navigation solutions. 
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Terrain models are also critical for mission planning of the airborne 
vehicles. Therefore, selection of the optimal path when using terrain models is one 
of the major problems of TAN. Paris and Le Cadre [24] investigate the planification 
of a mobile trajectory in order to use its own motion for improving its position 
estimation. In other words, an optimal trajectory is aimed to be planned which 
minimizes the localization error along the path or at the arrival area. In the paper of 
McFarland, et al [25], techniques originally developed for robot motion planning 
are applied to compute paths for autonomous air vehicles, such as cruise missiles or 
UAV’s. This approach is said to be particularly useful in multi-objective 
optimization problems such as intercepting a target while also maneuvering to 


minimize observability to ground-based tracking stations. 


Improvement of TAN using optimization is also one of the subjects of TAN. 
Bar-Gill, et al [26] propose a new method for improving the accuracy of TAN 
algorithms. They minimize the navigation errors which propagate along the flight 
path by designing airframe trajectories in a priori mission planning. The method 
uses information theory-based conditional entropy mapping and synthesizes 
minimum-entropy trajectories. Hence, by selecting optimal flight paths, navigation 


accuracies of the used TAN algorithms are improved. 


In the paper of Li, et al [27], optimal control methodology is adopted to 
design a terrain following controller for cruise missile. In this methodology, both 
tracking errors and control increments are considered in a guadratic penalty 
function. This paper is different from others; because, terrain following flight is 
investigated from the control point of view. Here, the TAN algorithms used and 


errors due to navigation system are not considered. 
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1.4.3. TAN Applications 


TAN applications can be found in literature especially for military purposes. 
Here, instead of developing new TAN algorithms, the applications of the known 
algorithms to real systems are given. In order to concentrate on TAN algorithms in 
detail, TAN applications are presented in a separate section. TAN algorithms found 


in literature will be investigated in the following section. 


Baird and Snyder [28] describe the design, mechanization and preliminary 
flight testing of a new altitude channel implementation, referenced primarily to the 
SITAN altitude estimates for AFTI/F-16 aircraft. Their paper is a typical SITAN 
algorithm application to a real system. In a similar way, Hollowell [29] presents the 
application of SITAN algorithm to US Army UH-1 Helicopter. In the paper, 
Multiple Model Adaptive Estimation (MMAF) techniques are employed for SITAN 
algorithm using a bank of single state Kalman filters to ensure that reliable position 


estimates are obtained even in the face of large initial position errors. 


Another example of a TAN application is the paper of Nordlund and 
Gustafsson [30]. They estimate the position of an aircraft using a terrain aided 
positioning algorithm based on a Rao-Blackwellisation technigue. This technigue 
uses recursive Monte Carlo methods, also known as particle filters and provides a 
favorable approximate solution. The TAN algorithm used here will be investigated 


in the following section considering the original paper of Bergman, et al [31]. 


TAN can be used in applications not only for air vehicles but for underwater 
and land vehicles as well. Newman and Durrant-Whyte [32] describe and 
investigate autonomous navigation of an underwater vehicle which uses inertial and 
sonar based sensors. In their paper, they associate inertial and sonar based world 
frame feature information in order to form a robust navigation algorithm. They do 
not use ready terrain map information for navigation; but, they form feature 


information around the underwater vehicle by using sonar measurements. 
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TAN applications for land vehicles are found in the papers of Madhavan, et 
al [33] and Bruder, et al [34]. Madhavan, et al [33] describe a TAN system which 
employs points of maximum curvature extracted from laser scan data as primary 
landmarks. On the other hand, Bruder, et al [34] present the development and 
implementation of a new sensor integration algorithm employing a terrain map to 
reduce INS errors. The algorithm used in this paper is very similar to classical TAN 


applications which use terrain height data for navigation correction. 


Terrain model improvement can also be considered as applications of TAN. 
Accurate terrain models not only improve the accuracy of the navigation system but 
they are also reguired for accurate height profile of the concerned areas. Morisue 
and Ikeda [35] demonstrate a navigation system which is used for high level of 
location accuracy. They achieve it by using various map-matching technigues. On 
the other hand, McLellan and Schleppe [36] describe an integrated real-time 
differential GPS and barometry system, with the prime aim of significantly 
changing and improving the method of positioning and layout of Shell Canada's 
land seismic surveys. They also state that the system had provided horizontal 
positioning better than 5 meters and height accuracy of better than 2 meters at 2 
sigmas. The system proposed is actually an integrated GPS system instead of a 
TAN algorithm application. However, since terrain height information is obtained 
accurately, the paper can be considered as an example for batch process terrain 


modeling application. 


1.4.4. TAN Algorithms 


As stated in the previous sections above, the hearth of TAN is the 
algorithms. A number of TAN technigues have been developed and tested. These 
fall into two general algorithmic categories of batch and recursive algorithms [11] 


as explained in the previous sections. 
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Literature survey on TAN algorithms is done considering these algorithmic 


categories. 


1.4.4.1. Batch TAN Algorithms 


The famous batch TAN algorithm found in literature is TERCOM. As it was 
mentioned before, original famous report of Baker and Clem (1977), named 
“Terrain Contour Matching (TERCOM) Primer” could not be obtained. However, 
detailed information about TERCOM is found in the book of Siouris [37]. 
TERCOM is a form of correlation guidance based on a comparison between the 
measured and the pre-stored features of the profile of the ground (i.e., terrain) over 
which a missile or aircraft is flying. Generally, terrain height forms the basis of this 
comparison [37]. There are a number of correlation algorithms (e.g., mean squared 
difference (MSD), mean absolute difference (MAD), the normalized MAD, the 
normalized MSD, and the product method) of varying complexity and accuracy that 
can be used to correlate the measured data with the reference data. Furthermore, the 
MAD algorithm provides the best combination of accuracy and computational 
efficiency for performing real-time terrain contour matching in an onboard 
computer environment [37]. Actually, TERCOM is a maximum likelihood 
estimator which uses only terrain height information for determining the vehicle's 
actual position. TERCOM is a batch process. Therefore, information about the 
position of the vehicle is post processed in order to have a navigation solution. 
TERCOM will be investigated in detail in the following chapter as one of the major 
TAN algorithms. 


Johnson, et al [11] improve the performance of TERCOM and SITAN by 
using maximum a posteriori estimator (MAP). Their technigue makes more 
complete use of the information content of the stored elevation data, the a priori 


knowledge of the errors in the stored elevation data and the elevation 
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measurements. Actually, the technigue is a batch algorithm which uses past 
information for TAN algorithm. However, this past information improves the 
results of both batch and recursive algorithms as explained in the paper. The theory 
behind the algorithm is straight forward; however, in order to apply MAP 


algorithm, extra computations are reguired. 


Erhui, et al [38] propose a new TAN algorithm based on the probability 
distribution differences of terrain height samples. They call their technigue as 
Probability-Based Terrain Aided Navigation (PTAN) approach. The technigue 
proposed is a batch algorithm and instead of correlating the height data collected by 
the radar and the barometer as in TERCOM, the proposed PTAN algorithm 
computes the probability distribution difference between them. The minimum 
probability distribution difference gives the best matching and the position of the air 
vehicle is determined accordingly. Again, the theory behind the algorithm is straight 
forward; however, in order to find probability distributions, considerable 


computational load is required. 


Zhou and Zhang [39] propose a scheme of TAN based on principle of 
computer vision. Being different from the conventional terrain matching technique, 
i.e. TERCOM, the scheme uses CCD camera rather than barometer and radio 
altimeter as sensing element. The technique proposed is a batch algorithm and it is 
claimed in the paper that shorter flight time is sufficient for successful terrain 
matching. Since, the original paper is in Chinese, details of the algorithm can not be 
obtained; only abstract of the paper is investigated. However, since CCD camera is 
used for correlation, environmental constraints should be considered. In other 
words, CCD camera can not be used in all weather conditions. This is thought to be 


the major drawback of the algorithm proposed. 


Quintang, et al [40] propose a new TAN approach using probabilistic data 
association filter (PDAF) to overcome irresolvable ambiguities in the correlation 


function used in TERCOM. The basic idea of the approach is to convert correlation 
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function value to the probability of position estimate being actual position of the air 
vehicle. It is shown via set of simulations that the method can improve the 
performance of TAN compared to TERCOM. The approach proposed is a batch 
algorithm. The interesting point is it uses one of the modem radar tracking 
algorithms. Actually, TAN is a data association problem, especially for the 
acguisition mode where INS position errors are very large. Here, it was thought 
whether the algorithm could be used for real-time applications. Therefore, this 
paper gave inspiration for implementing modern data association algorithms to 


TAN in the Ph D. study. 


1.4.4.2. Recursive TAN Algorithms 


The major recursive TAN algorithm found in literature is SITAN proposed 
by Hostetler and Andreas [12]. They investigate the application of nonlinear 
Kalman filtering techniques to the continuous updating of an INS using individual 
radar terrain clearance measurements in their paper. First order Extended Kalman 
Filter (EKF) is used in order to model the slopes of the terrain surface. Hence, real- 
time TAN solution can be obtained. Moreover, for large initial position 
uncertainties, a parallel Kalman filter technique which uses a bank of reduced order 
filters is used. The technique is the first EKF implementation to TAN. However, 
due to highly nonlinear structure of the terrain profiles, the filter solutions can 
diverge especially for large position errors. Moreover, linearization of the terrain 
slopes is the critical point for SITAN algorithm. Especially for mountainous 
terrains, modeling of the terrain slopes is a considerable problem. SITAN will be 


investigated in detail in the following chapter as one of the major TAN algorithms. 


Pei, et al [41] propose BITAN algorithm for navigation solution in their 
paper. BITAN is a type of TAN algorithm using the Kalman filtering theory to 


estimate position and velocity errors of the INS. The algorithm has been developed 
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for both acguisition and tracking modes of operation. Actually, the algorithm is an 
improved SITAN algorithm especially for acguisition mode. Moreover, again bank 
of Kalman filters are used for navigation solution. Hence, same problems for 


SITAN algorithm exist also for this algorithm. 


TAN is a nonlinear estimation problem. Bergman, et al [31] derive the 
optimal Bayesian solution for TAN. The implementation is grid based, calculating 
the probability of a set of points on an adaptively dense mesh. Actually, Bayes 
formula is a well-known formula in estimation. However, direct application of the 
formula is very restricted due to computational problems. As a result of this, 
Bergman, et al [31] propose the Cramer-Rao bound for Bayesian solution 
implementation. The major disadvantage of the algorithm proposed is originated 


from computational problems. 


One of the most interesting algorithms is proposed by Enns and Morrell 
[42]. They propose a new TAN algorithm called VATAN which uses the Viterbi 
algorithm for navigation solution. The Viterbi algorithm is a dynamic programming 
algorithm used for data association problem. From the simulation results, it is 
shown that VATAN algorithm overcomes divergence problems associated with the 
EKF in SITAN and provides position estimates with smaller average squared errors. 
Actually, navigation accuracy is improved with VATAN compared to SITAN 
especially for flat and mountainous terrains. In order to implement the algorithm, 
conditional probabilities of the measurements and INS states should be calculated 
recursively. Actually, the algorithm is different from other real-time TAN 


algorithms and it has better results than SITAN. 


Dezert [43] proposes a new application of PDAF for improving the accuracy 
of autonomous strapdown INS. The method proposed is a TAN algorithm based on 
landmark detection combined with a classical strapdown INS. It is also stated that 
the algorithm can be integrated with relatively low cost in existing operational TAN 


systems. Actually, the algorithm does not use elevation data for navigation solution. 
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However, it is a real-time application of PDAF and relation with the former paper 
of Ouintang, et al [40] can be obtained where batch implementation of PDAF is 
used. Therefore, this paper also gave inspiration for implementing real-time PDAF 


to TAN in the Ph.D. study. 


Some recursive TAN algorithms are also proposed which use images for 
navigation solution. Hongbo, et al [44] and Bevington, et al [45] can be examples to 
image based TAN. Hongbo, et al [44] propose a TAN algorithm which use range 
images from imaging laser radar. On the other hand, Bevington, et al [45] use 
images of Synthetic Aperture Radar (SAR) for navigation solution. Both methods 
reguire detection of land marks since images are used. As it was stated in the 
previous section, the major drawback of the image based TAN comes from 


environmental constraints. 


There exist also some hybrid TAN methods which use both batch and 
recursive algorithms together. Metzger, et al [3] propose a hybrid TAN system 
which uses a bank of Kalman filters and a comparison technigue. Actually, the 
proposed algorithm is a mixture of TERCOM and SITAN algorithms. Using the 


advantages of both algorithms, better navigation solutions can be obtained. 


1.4.5. TAN Patents 


United States patents related with TAN are also investigated in literature 
survey. Since patents are practical applications, detailed information can be 


obtained from them. Several patents are investigated related with TAN algorithms. 


Chan and Snyder [46] propose a system for correlation and recognition of 


terrain elevation. They use correlation function in freguency domain in order to 
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improve navigation solutions. The method is the improvement of the correlation 


technigue. 


Baird [47] integrates TERCOM and SITAN algorithms with a modified 
Kalman filter processor. Hence, the operation of the SITAN processing is 
effectively continuously optimized. Actually, the system is the application of a 


hybrid TAN algorithm as explained in the previous section. 


Lerche [48] improves TERCOM method by scanning a larger area for 
correlation process. The method is applied for the navigation of an aircraft. Raymer, 
et al [49] proposes a method for Schuler cycle error reduction for use in a TAN 


system. By detecting Schuler cycles, TAN system errors are degraded. 


Finally, Goebel, et al [50] propose a terrain correlation system for TAN. 
Actually, the correlation system is an improved TERCOM algorithm which uses 
MAD correlation. Detailed information including application methods are given in 


the related patent. 


It is known that TAN algorithms are used generally for military purposes. 
Due to the confidentiality of the subject, related patents about TAN are limited. In 


fact, the patents found are taken many years later than their technology developed. 


1.5. Target Tracking 


TAN is a nonlinear estimation problem; since, terrain height information is 
used for navigation solution. Actually, TAN can be considered as a data association 
problem, especially for the acguisition operation mode where INS position errors 
are considerably large. From the literature survey of Ouintang, et al [40] and Dezert 


[43], it has been thought that modern data association algorithms can be 
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implemented for real-time TAN algorithms. Therefore, radar tracking, especially 
data association subject is investigated. In this section, general information about 


radar tracking and data association algorithms will be given. 


1.5.1. Background 


The modern need for tracking algorithms began with the development of 
radar during World War II. By the 1950’s, radar was a relatively mature technology. 
Systems were installed aboard military ships and aircraft and at airports. The 
tracking of radar targets, however, was still performed manually by drawing lines 
through blips on a display screen. The first attempts to automate the tracking 
process were modeled closely on human performance. For the single-target case, 
the resulting algorithm was straight forward; the computer accumulated a series of 
positions from radar reports and estimated the velocity of the target to predict its 


future position (51). 


Even single-target tracking presented certain challenges related to the 
uncertainty inherent in position measurements. A first problem involves deciding 
how to represent this uncertainty. A crude approach is to define an error radius 
surrounding the position estimate. This practice implies that the probability of 
finding the target is uniformly distributed throughout the volume of a three- 
dimensional sphere. Unfortunately, this simple approach is far from optimal. The 
error region associated with many sensors is highly non-spherical; radar, for 
example, tends to provide accurate range information but has relatively poorer 
radial resolution. Furthermore, one would expect the actual position of the target to 
be closer on average to the mean position estimate than to the perimeter of the error 
volume, which suggests, in turn, that the probability density should be greater near 


the center [51]. 
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A second difficulty in handling uncertainty is determining how to interpolate 
the actual trajectory of the target from multiple measurements, each with its own 
error allowance. For targets known to have constant velocity (e.g., they travel in a 
straight line at constant speed), there are methods for calculating tile straight-line 
path that best fits, by some measure, the series of past positions. A desirable 
property of this approach is that it should always converge on the correct path, as 
the number of reports increases, the difference between the estimated velocity and 
the actual velocity should approach zero. On the other hand, retaining all past 
reports of a target and recalculating the entire trajectory every time a new report 
arrives is impractical. Such a method would eventually exceed all constraints on 


computation time and storage space [51]. 


A near-optimal method for addressing a large class of tracking problems was 
developed in 1960 by R.F. Kalman. His approach, referred to as Kalman filtering, 
involves the recursive fusion of noisy measurements to produce an accurate 
estimate of the state of a system of interest. A key feature of the Kalman filter is its 
representation of state estimates in terms of mean vectors and error covariance 
matrices, where a covariance matrix provides an estimate (usually a conservative 
over-estimate) of the second moment of the error distribution associated with the 
mean estimate. The sguare root of the estimated covariance gives an estimate of the 
standard deviation. If the seguences of measurement errors are statistically 
independent, the Kalman filter produces a seguence of conservative fused estimates 


with diminishing error covariances [51]. 


Kalman’s work had a dramatic impact on the field of target tracking in 
particular and data fusion in general. By the mid-1960’s, Kalman filtering was a 
standard methodology. It has become as central to multiple-target tracking as it has 
been to single-target tracking; however, it addresses only one aspect of the overall 


problem [51]. 
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1.5.2. Data Association Algorithms 


In tracking targets with less-than-unity probability of detection in the 
presence of false alarms (clutter), data association, deciding which of the received 
multiple measurements to use to update each track is crucial. A number of 


algorithms have been developed to solve this problem. Two simple solutions are; 


1. Strongest Neighbor Filter (SNF), and 


2. Nearest Neighbor Filter (NNF). 


In the SNF, the signal with the highest intensity among the validated 
measurements (in a gate) is used for track update and the others are discarded. In 
the NNF, the measurement closest to the predicted measurement is used. While 
these simple technigues work reasonably well with benign targets in sparse 
scenarios, they begin to fail as the false alarm rate increases or with low observable 


(low probability of target detection) maneuvering targets [52]. 


The NNF is perhaps the simplest approach for determining which tracked 
object produced a given sensor report. When a new position report arrives, all 
existing tracks are projected forward to the time of the new measurement. Then, the 
distance from the report to each projected position is calculated, and the report is 
associated with the nearest track. More generally, the distance calculation is 
computed to reflect the relative uncertainties (covariances) associated with each 


track and report [51]. 


In Figure 8, NNF implementation is shown. The idea of the rule is to 
estimate each objects position at the time of a new position report, and then assign 
the report to the nearest such estimate. This intuitively plausible approach is 
especially attractive because it decomposes the multiple-target tracking problem 


into a set of single-target problems [51]. 
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Figure 8. NNF Implementation (51) 


Data association becomes more difficult with multiple targets where the 
tracks compete for measurements. Here, in addition to a track validating multiple 
measurements as in the single target case, a measurement itself can be validated by 
multiple tracks (i.e., contention occurs among tracks for measurements). Several 


algorithms are developed to handle this contention: 


1. Track Splitting (TS), 


2. Multiple Hypothesis Tracking (MHT), 


3. Probabilistic Data Association (PDA), 


4. Joint Probabilistic Data Association (JPDA). 
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Actually, there are various data association algorithms proposed for multiple 
target tracking in literature. However, the algorithms listed above can be considered 


as modern data association algorithms. 


In track splitting and MHT, a robust solution to the problem of assignment 
ambiguities is found by creating multiple hypothesis tracks. Under this scheme, the 
tracking system does not have to commit immediately or irrevocably to a single 
assignment of each report. If a report is highly correlated with more than one track, 
an updated copy of each track can be created; subseguent reports can be used to 
determine which assignment is correct. As more reports come in, the track 
associated with the correct assignment will rapidly converge on the true target 
trajectory, whereas the falsely updated tracks are less likely to be correlated with 


subsequent reports [51]. 


This basic technique is called track splitting. One of its worrisome 
consequences is a proliferation in the number of tracks upon which a program must 
keep tabs. The proliferation can be controlled with the same track deletion 
mechanism used in the nearest-neighbor algorithm, which scans through all the 
tracks from time to time and eliminates those that have a low probability of 
association with recent reports. A more sophisticated approach to track splitting, 
called multiple-hypothesis tracking, maintains a history of track branchings, so that 
as soon as one branch is confirmed, the alternative branches can be pruned away 
[51]. MHT is a more powerful (but much more complex) algorithm that handles the 
multi-target tracking problem by evaluating the likelihood that there is a target 


given a sequence of measurements. 


In PDA, instead of using only one measurement among the received ones 
and discarding the others, an all of the validated measurements with different 
weights (probabilities) are used. The standard PDA and its numerous improved 
versions have been shown to be very effective in tracking a single target in clutter 


[52]. 


a 


JPDA algorithm is used to track multiple targets by evaluating the 
measurement-to-track association probabilities and combining them to find the state 
estimate [52]. Actually, JPDA is the developed version of PDA algorithm for 


multiple targets. 


PDA, TS and MHT will be investigated in detail in the following chapters 


for TAN implementation. 


1.6. Outline of the Thesis 


In this section, the outline of the Ph.D. study will be given. The thesis is 


composed of five chapters: 


1. Introduction, 


2. Major TAN Methods, 


3. Implementation of Radar Tracking Algorithms to TAN, 


4. Case Study, 


5. Discussion and Conclusion. 


In the first chapter, an introduction to the study was done. First, the scope of 
the study was presented. Then, general information about cruise missiles and TAN 
was given. A detailed literature survey was performed about TAN and was 
presented in this chapter. Finally, general information about data association 


algorithms was given as fundamental knowledge of the study. 
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In Chapter 2, major TAN methods are investigated. First, INS errors of the 
cruise missiles and need for TAN systems are discussed. Then, major TAN methods 
including TERCOM, SITAN and VATAN are presented in detail. Fundamentals of 
the major methods are discussed in this chapter in order to make comparisons for 


the implemented TAN algorithms in the Ph.D. study. 


In Chapter 3, implementation of data association algorithms to TAN is 
presented. This chapter contains the original Ph.D. work. First, general information 
about modem target tracking algorithms are given. PDAF and TSF algorithms and 
their general implementations are investigated. Then, PDAF and TSF 
implementations to TAN are presented. At the end of the chapter, a simple 
simulation model is developed for the mid-course flight of the cruise missile. 
Finally, simulations are performed with the implemented TAN algorithms and the 


results are compared with the major TAN methods. 


In Chapter 4, case studies are performed. A 6 DOF simulation tool is 
developed for the simulation of the mid-course flight of a cruise missile. 
Implemented TAN algorithms are used with the 6 DOF simulation model and their 


performances are investigated. 


In Chapter 5, the results obtained from the study are discussed. Advantages 
and disadvantages of the new implemented TAN algorithms are compared with the 


major TAN algorithms. Finally, conclusions of the study are presented. 
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CHAPTER 2 


MAJOR TERRAIN AIDED NAVIGATION METHODS 


In this chapter, first, general INS errors and TAN INS errors in cruise 
missiles will be discussed. Then, the need for TAN in cruise missiles will be 
investigated. Next, major TAN algorithms and their implementations will be 
presented. TERCOM, SITAN and VATAN will be investigated in detail in this 
chapter. Eventually, navigation performance of these major TAN algorithms will be 


discussed; and, conclusions obtained will be presented. 


2.1. Cruise Missile INS Errors 


2.1.1. INS Only Errors 


The development of inertial navigation technology took place primarily in 
Germany, the United States and the former Soviet Union. The gyro compass 
indicating true north on a moving base as on ships can be regarded as the beginning 
of inertial navigation. At the end of World War I the allies had in the Treaty of 
Versailles imposed restrictions to Germany for the maximum size of ships to be 
built. These restrictions promoted in this country gun stabilization and inertial 


technology in general, which culminated at the end of World War II in a 


40 


functioning air-supported gyrocompass with electronic Schuler tuning for the “One- 
Man Submarines”, in the V2 guidance system and a true concept for an INS. After 
the war the development of this technology was taken over by the superpowers, the 


United States and the former Soviet Union [53]. 


Inertial navigation systems (INS) are sophisticated autonomous, 
electromechanical systems that supply the position, velocity and attitude of the 
vehicle on which they are mounted. INS is basically a measuring system; therefore, 
the outputs of an INS will contain errors due to its sensors (accelerometers and 
gyroscopes) and mechanization. Inertial navigation sensor component errors create 
error in the navigation system’s computed position, velocity, and attitude. 
Accelerometer and gyroscope errors can be represented in a general form, including 


some significant environment dependent errors as [54]: 


1. Biases and drifts, 


2. Scale factor and misalignments, 


3. White noise, 


4. Time correlated short-term errors, 


5. Other environment sensitive errors. 


A lot of experience has been gained on the behavior of INS errors from the 
accumulated experience of INS users and analysts. Various linear models were 
developed that describe accurately the behavior of these errors as given in 
references [54], [55], and [56]. These models were used in the implementation of 


Kalman filters for estimating the INS error outputs and error sources [55]. 
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There are two approaches to the derivation of INS error models. One of 


them is known as the ¢-angle (perturbation or true frame) approach, and the other 
is known as w-angle (or computer frame) approach. When deriving the 


perturbation error model, the nominal non-linear navigation eguations are perturbed 
in the local-level north-pointing Cartesian coordinate system that corresponds to the 


true geographic location of the INS. The y -angle error model, on the other hand, is 


obtained when the nominal equations are perturbed in the local-level north-pointing 
coordinate system that corresponds to the geographic location indicated by the INS. 
It has been shown that both models are equivalent and yield, therefore, identical 
results. The differential equations that describe the error behavior of the INS are 
divided into equations describing the propagation of the attitude errors. Both the 
translatory and the attitude error equations can be expressed in two different ways 
that yield two versions of the translatory error equations and two versions of the 
attitude error equations. The two versions of the translatory equations depend on 
whether the equation variables are position error components or velocity error 
components. The two versions of the attitude equations depend on whether the 
equation variables are components of the platform to computer frame attitude 
difference, or components of the platform to true frame attitude difference. All these 
versions are, of course, identical. In order to obtain a complete set of INS error 
equations, the analyst has to decide whether to adopt the perturbation or wy -angle 
approach. Once this choice is made, the analyst has to decide which of the two 
corresponding versions of the translatory equations to use and which of the two 


versions of the attitude equations to use. (These two choices are independent.) [56]. 


Most of the published work on INS errors adopt they -angle approach and 
use the velocity error version of the translatory error equation. This model is also 
used in the present analysis. In addition, the components of the platform to 
computer frame attitude differences are used as the variables of the attitude error 
equations. Although this angular difference is imaginary and cannot be measured, it 


possesses the advantage that the translatory error is not coupled into the attitude 
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error eguations. The physical attitude difference between the platform and the local- 
level north-pointing coordinate system is calculable using the position and attitude 
errors obtained from the solution of these INS error eguations. Then, a complete 


terrestrial INS error model, expressed by the following eguations is obtained [55]: 


OV +(0+0)x 07 =V-0Pxf+45 


(2.1) 
OF + px OF = dv (2.2) 
OV +x dP =e (2.3) 


where dv, 67 and OV are, respectively, the velocity, position, and attitude 
error vectors; Q is the Earth rate vector; © is the angular rate vector of the true 
coordinate system with respect to inertial frame; V is the accelerometer error 
vector; f is the specific force (accelerometer readings) vector; Ag is the error in 
the computed gravity vector; Ø is the vector of the rate of turn of the true frame 


with respect to Earth; and finally € is the gyro drift vector. From geometric 
relations, it can be shown that in the local north, east and down coordinate system 


(i.e. in the geographic frame) [55]: 


O cos 
0 (2.4) 
-Q sin A 


OI 
I 


where 4 is the local latitude. In the same manner, the vector © is computed 


as follows [55]: 
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ö-04p (2.5) 


fi-cosa 
where p = SÄ and u is the local longitude. 
—f4-sind 


When INS position, velocity and attitude error equations are resolved in true 


frame (i.e. geographic frame), nine scalar differential equations are obtained, which 


can be put in a state-space model. If the expressions for Q, & and 5 are used, the 


resulting state-space model is obtained as follows [55]: 
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(2.6) 
where; 
ör: Scalar position errors 


öv: Scalar velocity errors 


ow: Scalar attitude errors 


©:  Earth’s inertial angular velocity 
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H: Longitude of the true frame with respect to Greenwich meridian 


HK: Longitude rate of the true frame 
A: Latitude of the true frame with respect to Eguator 
À: Latitude rate of the true frame 


g: Farth's gravity 


R: Radius of Farth 


S: Specific forces sensed by the accelerometers 


V: Scalar accelerometer biases 


E: Scalar gyro drifts 


N, E, D : Subscripts denoting north, east and down components respectively 


s: Sine of the defined angle 


a Cosine of the defined angle 


The error model given in eguation (2.6) can be used in simulations for 
predicting INS errors of the system. Another way of determining INS errors is the 
direct application of the real error sources (from both sensors and mechanization) in 


the navigation eguations. Actually, the error model obtained above is used for 
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integrating various navigation systems with INS. However, for simulations, INS 


error model will be sufficient and its implementation will be much easier. 


INS error model is applied for the mid-course phase of a cruise missile. The 
cruise missile is assumed to be moving with constant velocity. Moreover, 
acceleration changes during mid-course flight are assumed smaller. INS guality is 


taken as 1.0 nm/hr for simulations. 


INS guality is the major parameter in order to achieve the reguired 
navigation solutions. It is mainly determined by sensor guality and initial alignment 
errors. In Table 5, various INS gualities and corresponding sensor and initial 


alignment errors are presented. 


Table 5. INS Sensor Error Sources [57] 


















































INS Quality (All errors except random walk are lo biases) 
Error Source 10 nm/hr 1.0 nm/hr 0.5 nm/hr 0.2 nm/hr 
Accelerometer Bias |223 ug 37 ug 19 ug 4.2 ug 
Accel. Scale Factor | 223 ppm 179 ppm 90 ppm 21 ppm 
Input Axis Misalign. | 22 arcsec 3 arcsec 1.5 arcsec 0.4 arcsec 
Random Walk 56 ug Nhz 56 ug Nhz 7.5 ug/Nhz 4.2 ug/Vhz 
Gyro Bias 0.11 deg/hr 4.5e-3 deg/hr |2.2e-3 deg/hr | 8.4e-4 deg/hr 
Gyro Scale Factor 112 ppm 112 ppm 7.5 ppm 1.67 ppm 
Input Axis Misalign. | 22 arcsec 2.2 arcsec 1.1 arcsec 0.4 arcsec 
Random Walk 0.078 deg/Vhr | 2.2e-3 deg/Vhr | 1.1e-3 deg/Vhr | 5e-4 deg/Vhr 
Initial Misalignment | 2089 arcsec/ o | 606 arcsec/ 600 arcsec/ 600 arcsec/ 
(Vertical/Horizontal) | 59 arcsec 59 arcsec 29 arcsec 29 arcsec 
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INS guality is generally expressed by the total position error divided by 
time. For example, 10 nm/hr INS guality will be sufficient for ballistic missiles. 
However, for military aircrafts and cruise missiles 1.0 nm/hr INS guality is 
reguired. In the same manner, as operation time and reguired range increases, INS 
guality will be also increased. Intercontinental Ballistic Missiles (ICBM) and space 
vehicles use very accurate INS. The major problem of using very accurate INS is its 
cost. Moreover, due to large space reguirements of very accurate INS, they can not 


be used in most of the military systems. 


INS error model simulations are performed in Simulink [58]. Considering 


1.0 nm/hr INS guality, horizontal position and velocity errors and attitude errors are 
obtained. It is known that, an initial altitude error (Ah, ) or altitude-rate error (Ah,) 


or an accelerometer error will grow exponentially with time, thus making the 
indicated altitude and altitude-rate indications useless after a few minutes. The 
instability of the vertical channel for INS will result, no matter how carefully the 
vertical component of gravity is mechanized as a function of computed altitude 
[59]. In real systems, vertical channel of INS is generally aided by barometric 
altimeters. Since, horizontal position errors are critical for TAN, altitude errors are 


not investigated. Simulation results are shown in Figure 9, Figure 10, and Figure 11. 
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Horizontal INS Position Errors vs. Time 
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Figure 9. Horizontal Position Errors of the INS Error Model 


Horizontal INS Velocity Errors vs. Time 
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Figure 10. Horizontal Velocity Errors of the INS Error Model 
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x 105 INS Attitude Errors vs. Time 
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Figure 11. Attitude Errors of the INS Error Model 


As it can be seen from the simulation results, guadratic increase of position 
errors due to double integration dominates the INS only navigation solution. On the 
other hand, the validity of the error model can be seen from Figure 9. Here, using 


sensor errors defined in Table 5, 1.0 nm/hr guality INS is achieved. 


2.12. TANINS Errors 


INS only errors of a navigation system are discussed in the previous section. 
As it can be seen from the simulation results, due to large navigation times of cruise 
missiles, INS should be aided with other navigation systems. TAN is the well- 


known method for improving navigation solution. 


The accuracy of the TAN position estimate for a simple case is derived from 


application of linear estimation theory. Using horizontal INS position errors 
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modeled as independent random walks (uncorrelated white noises), following 
expression can be obtained for the circular error probable (CEP) of horizontal 


position updates [13]. 


CEP, =0.57-6V"*.(Ad/s)”* -(o,/h)"* (2.7) 


where; 


CEP, : The steady state CEP of horizontal position updates (m), 


SS 


0, : Standard deviation of the profile measurement errors (m), 


h: Deterministic local terrain slope at the measurement locations in both 


down-range and cross-range directions (unitless), 
Ad: Distance between profile measurements (m), 
s: Vehicle ground speed (m/s), 


öV: Maximum INS velocity error (m/s). 


The primary value of the equation above is that it shows the sensitivities of 
accuracy to implementation parameters. Steady state CEP is most sensitive to the 


ratio o/h, least sensitive to ÖV , and nominally sensitive to the time between 


profile measurements Ad /s . Using typical values of; 
OV = 1 m/s (1 nm/hr-class INS) 
s =250 m/s 
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Ad = 100m 


h =0.05 (moderately rough terrain) 


o =15m 


results ina CEP of 29 m. Because of the assumptions in the equation given 


above, predictions should be treated as approximations, a conservative lower bound 


for TAN accuracy [13]. 


As it can be seen from linear TAN estimation results, INS error growth in 
time is limited using TAN algorithms. By correlating terrain profiles with INS 
solutions a few times during operation or recursively, position estimates are 


obtained. Then, INS is updated according to the estimated navigation solutions. 


22. TERCOM 


2.2.1. Background 


Terrain Contour Matching (TERCOM) can be defined as a technique for 
determination of the position location of an airborne vehicle with respect to the 
terrain over which the vehicle is flying. More specifically, TERCOM is a form of 
correlation guidance based on a comparison between the measured and the pre- 
stored features of the profile of the ground (i.e., terrain) over which a missile or 
aircraft is flying. Generally, terrain height forms the basis of this comparison. 
Reference terrain elevation source data descriptive of the relative elevations of the 


terrain in the fix point areas are stored in the air vehicle’s onboard computer. 
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Obtaining the reference data reguires prior measurement of the ground contours of 
interest. These data are in the form of a horizontally arranged matrix of digital 
elevation numbers. A given set of these numbers describes a terrain profile. The 
length of contour profile necessary for a unigue fit is a function of terrain 


roughness, but is in the range of 6 to 10 km and can be a curved path [37]. 


As the vehicle flies over the matrix area, data describing the actual terrain 
profile beneath the vehicle are acguired. That is, the actual profile is acguired using 
a combination of radar and barometric altimeter outputs sampled at specific 
intervals, and when compared against the stored matrix profiles provide the position 
location. This type of guidance is used for updating a mid-course guidance system 
on a periodic basis, and has been applied to the guidance of cruise missiles, which 
usually fly at subsonic speeds and fairly constant altitude. With regard to mid- 
course guidance, it is well known that the simplest mid-course guidance is the 
explicit guidance method. The guidance algorithm has the capability to guide the 
missile to a desired point in the air while controlling the approach angle and 
minimizing an appropriate cost function. Furthermore, the guidance gains of the 
explicit guidance law are usually selected to shape the trajectory for the desired 


conditions [37]. 


The TERCOM technigue, first patented in 1958, relies for its operating 
principle on the simple fact that the altitude of the ground above sea level varies as 
a function of location. Historically, TERCOM has evolved from several R&D 
programs that developed certain areas of the overall process. These programs 
perfected the technology as it is known today [37]. In Table 6, a chronological 


overview of this development is summarized. 


In the following sections, TERCOM method will be investigated in detail. 
Fundamentals of the TERCOM concept are taken from Siouris [37]. At the end of 


TERCOM section, simulations performed will be discussed. 
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Table 6. Chronological Overview of TERCOM Development [37] 












































PROGRAM YEAR OBJECTIVES 
Guidance package for SLAM 
missile TERCOM concept 

Fingerprint 1958 first 
proposed. 
Feasibility study of terrain 
TERCOM 1960-1961 contour matching. 
Design and development of a 
LACOM (Low Altitude 1963-1965 complete fix-taking 
Contour Matching) subsystem. 
Improve TERCOM 
RACOM (Rapid Contour computation procedures and 
: 1963-1966 : 
Matching) increase accuracy. 
SAMSO (USAF’s Space 
and Missiles Systems 
Organization) Programs 
(a) TPLS (Terminal 
Position Location System Application of terrain 
(b) TERSE (Terminal 1963-1971 correlation technigues for 
Sensing Experiment) N ballistic missiles. 
(c) TERF (Terminal Fix). 
(d) TSOFT (Terminal 
Sensor Overland Flight 
Test). 
Study and define a TERCOMI 
poni drone system capable of 
Avionics Update 1972-1975 operational deployment. 
Feasibility study for 
POIS incorporation in cruise missile 
OT (Terrain Aided 1979-1974 and evaluation of snow 
TW 972-197 coverage effects on terrain 
profile acquisition. 
McDonnell — Douglas 
Astrodynamics awarded a 
Competitive Flyoff 1975 contract for TERCOM 
system. 
RACOM (Recursive All Improve terrain correlation 
Weather Contour 1975 update accuracy. 
Matching) 
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2.2.2. TERCOM Concept 


TERCOM system uses an airborne altimeter and a data processor to 
correlate the measured terrain contours to obtain the best estimate of position. The 
TERCOM system relies on a set of digital maps stored in the memory of the 
missile’s onboard computer. These maps consist of rectangular arrays of numbered 
squares representing the variation of ground elevation above sea level as a function 
of location. Consequently, as the missile approaches an area for which the computer 
memory has a map, the onboard radar altimeter starts providing a stream of ground 
elevation data. Furthermore, the computer, by comparing these data with the 
information it has in its memory, can accurately determine the actual trajectory of 


the missile and instruct the autopilot to return the missile to its planned trajectory. 


Four such corrective maneuvers are shown in the vertical overhead view in Figure 


12 [37]. 
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Figure 12. TERCOM Maps in Use [37] 


The map types used in TERCOM differ in length, width, and cell size. The 
cell size determines, in part, the accuracy of the TERCOM fix. The TERCOM maps 
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become smaller and are spaced closer together as the missile approaches the target. 
As a result, because of the decreasing cell size, the updates become more accurate. 
A terminal accuracy on the order of 100 meters (i.e. DTED Level 1) is considered 


feasible for the TERCOM system [37]. 


The process of determining air vehicle position by the use of terrain contour 
matching can generally be described as consisting of three basic steps; data 
preparation, data acguisition, and data correlation. In Figure 13, TERCOM concept 


is illustrated [37]. 
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Figure 13. TERCOM Concept [37] 


The critical part of TERCOM process is the data correlation where 
navigation solution is performed. There are several data correlation algorithms like 


MAD and MSD and they will be discussed in the following section. 
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TERCOM measurement process is illustrated in Figure 14 in block diagram 
form. The radar altimeter acguires altitude estimates above terrain. Then, the radar 
altimeter output is differenced with the system's reference altitude. Various 
arithmetic operations (e.g. mean removal and guantization) are then performed on 
the differenced data. Finally, the correlation between the stored and acguired data is 


performed with the MAD function, and a position fix is determined (37). 


Terrain 
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Figure 14. TERCOM Measurements [37] 


2.2.3. TERCOM Data Correlation Technigues 


There are a number of correlation algorithms (e.g., mean sguared difference 


(MSD), mean absolute difference (MAD), the normalized MAD, the normalized 
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MSD, and the product method) of varying complexity used in TERCOM. 
Furthermore, the MAD algorithm provides the best combination of accuracy and 
computational efficiency for performing real-time terrain contour matching in an 
onboard computer environment. Therefore, here only the MAD and MSD 


correlation algorithms will be discussed [37]. 


The MAD algorithm is applied considering the first N height differences to 
be acguired. Then, these differences are removed, so that the sample profile is its 
mean value. Next, this profile is compared with each row of matrix data in the 
following manner. Let h, (1<n<N) denote any row of matrix data and H, the 
seguence of reguired data. Conseguently, the MAD algorithm, which is used for 
correlating the measured terrain elevation file with each down-track column of the 


reference matrix, is defined as follows [37]: 


MAD, ,, = (Yh n -H 
izl 





(2.8) 


m,n 


where; 

MAD, ,, : The value of the mean absolute difference between the k’th terrain 
elevation file and the m’th reference matrix column, 

N: The number of samples in the measured terrain elevation file and 
usually it is also equal to the number of rows in the reference 
matrix, 

M: The number of reference matrix columns, 

K: The number of measured terrain elevation files used in the 


correlation process, 
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| |= The absolute value of the argument, 


n,m,k: Row, column, and terrain elevation file indices, 





H : The stored reference matrix data, 1<m <M,1<n<N, 


The k’th measured terrain elevation file, 1<k <K. 


The MSD algorithm can be expressed in terms of the profile in question. 


Mathematically, the expression for MSD is [37], 
= 2 
MSD, =(1/ MY (S; -5;) (2.9) 
i=1 


where, 


S,,S,:J’th and k’th profiles, 


N: Length of each profile. 


Note that for uniformity, the MAD algorithm can also be expressed as in the 
expression for the MSD. Thus, 


N 
MAD, = VM), 





Si — Six (2.10) 
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Examination of the expressions for the MAD and MSD processors indicates 
that both of these correlators can be viewed as distance measures, where the 
dimensions of the space for which these distances are defined correspond to the 
number of elements in the profiles. From (2.9) and (2.10), it is noted that the 
ambiguity between any two profiles is defined as the probability “ P ” that sensed 
data corresponding to one of the profiles will be closer (in terms of the distance 


measure) to the other profile than to the one from which it was taken [37]. 


Mathematically, the ambiguity € can be expressed as: 


E a jk < Cy], Where a minimum of C, is sought, | 
ik = 


PIC,, >C,,], where a maximum of C., is sought. (2.11) 
For a MAD processor, C is given by the following expression: 
N 
rr SS (2.12) 
i=1 





where, 


S,: o J'th measured profile, 


R,: kth reference profile. 


A more detailed account of the terrain correlation processing for a single map 


is conceptually shown in Figure 15. 
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Figure 15. Terrain Correlation Processing [37] 


2.2.4. Terrain Roughness Characteristics 


For TERCOM correlation process, roughness and unigueness of the selected 
terrain is very critical. It should be noted that the TERCOM concept will not work 
over all types of terrain. For instance, the rougher the terrain, the better TERCOM 
works. However, good terrain must be more than just rough, it must be unique (i.e., 


a given profile out of the TERCOM map must not resemble any other map [37]. 


Terrain roughness is defined as the standard deviation of the terrain 
elevation samples as shown in Figure 16. It is usually referred to as “sigma-T” (or 


07) [37]. 


Sigma- T is defined by the eguation: 
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O= lam H, -HY (2.13) 


where, 


N 
H =(1/N i H,: Mean Elevation 


izl 


Mean elevation Terrain profile 
| i 





Mean sea level 


Figure 16. Terrain Standard Deviation (Sigma-T) [37] 


Thus, o, is a measure of the variation of the terrain elevation about its 
average elevation. Note that the minimum value of o, required to support 
TERCOM operation is approximately 25 ft (7.62 m). Areas that have sigma-T 
values of fifty or greater are usually considered as good candidates for TERCOM 
fix areas. Obviously, lakes and very flat or smooth areas have low values of sigma- 
T. Therefore, they are not suitable as fix areas. However, sigma-T is not the only 


criterion for determining whether a given area is suitable for TERCOM operation 


(371. 


In particular, there are three parameters that are used to describe TERCOM- 
related terrain, and their values can give an indication of the terrain's ability to 


support a successful TERCOM fix. These parameters are sigma-T, sigma-Z (o,), 
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and the terrain correlation length ( X,.). It is usually assumed that parallel terrain 


elevation profiles that are separated by a distance greater than X, are independent 


of each other [37]. 


Sigma-Z is defined as the standard deviation of the point-to-point changes in 
terrain elevation (i.e., the slope) as shown in Figure 17. Like sigma-T, the value of 
sigma-Z provides a direct indication of terrain roughness. Sigma-Z has also been 
shown to be a valid indicator of TERCOM performance. The expression for sigma- 
Z, assuming a Gaussian autocorrelation function, can be obtained from Figure 17. 


Mathematically, sigma-Z is given by the eguation [37]: 





o, = Jou — DIY (D,-D)? (2.14) 
where, 
D=H-H 


D=(1I/(N - Dy D, 


i=1 


Terrain profile 
/ 






Mean sea level 


Figure 17. Definition of Sigma-Z [37] 
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The two parameters sigma-T and sigma-Z are related to the third parameter 


X, according to the relation [37]: 


o, =2-07" -[1-exp(-Ad/X,)'] (2.15) 
where; 
Ad: Cell size (or distance between elevation samples). 


2.2.5. Simulations and Discussion 


In order to investigate TERCOM performance, a simulation model is 
developed with Matlab [60]. The sample map considered for the simulations has the 
size of M=21 by N=100 where the cells are 100 x 100 meters approximately. The 
TERCOM procedure is as follows: 


1. Map is selected considering CEP of the INS. 


2. MxM (21 x21) for N=1 is considered. 


3. Height measurements are considered then. 


4. Absolute differences, n= 1 to M and m= 1 to M are 





h,, <H 





m,n 


calculated for MAD process. (21 x 21 operations) 


—H, ) n= 1 to M and m= 1 to M are 


1,m m,n 


5. Square differences, (A 


calculated for MSD process. (21 x 21 operations) 
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6. Steps “1” to “3” are repeated for N=1 to 100. (21 x 21 x 100 


operations) 


is calculated for MAD process. 


m,n |? 


7. MAD,,= UNY Vin -H 
izl 


(Extra sum and averaging operations for 21 x 21 x 100 elements) 


N 
8. MSD, = (/N)Y (S, =S y , 1s calculated for MSD process. (Extra 
izl 


sum and averaging operations for 21 x21 x 100 elements) 


9. Minimum of MAD and MSD functions are sought in order to 
determine the indices “i" and “j” of the horizontal position fixes for 
both MAD and MSD processes. (Determination of the minimum 
points) 


For the simulations, DTED Level 1 data were reguired and they have been 
obtained from HGK. The properties of DTED prepared for Turkey were given in 
Table 3 [16]. Horizontal accuracy of Level 1 DTED is defined as +130 m, and 
vertical accuracy as +30 m. In fact, especially horizontal accuracy of the DTED 
Level 1 data for Turkey is not sufficient for navigation purposes. Actually, DTED 
Level 2 data which have horizontal accuracy of +26 m can be resampled to DTED 
Level 1 and used for practical applications. However, for the Ph.D. study, DTED 


Level 1 data are used considering the horizontal accuracies of DTED Level 2. 


In order to perform the simulations, first selection of the areas is performed 
using the mapping software OziExplorer [61]. In order to select the areas for 
simulation, roughness of the surfaces is investigated using elevation property of the 
software which uses DTED Level 1 files obtained from HGK. Sample area 


selection using OziExplorer is shown in Figure 18. 
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Figure 18. Sample Area Selection from OziExplorer Software [61] 


Then, using OziExplorer3D [62] software, optional add-on to the 
OziExplorer software which allows map images to be viewed in 3D, selected areas 
were rendered as matrix grids. For the simulations, three special areas were 


selected: 
1. Area with rough surface, 
2. Area with smooth surface, 


3. Area with having uniqueness (i.e. a single mountain). 
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, M x N area 


Selected areas rendered using OziExplorer3D are shown in Figure 19 and 


Figure 20. 






























Figure 19. Rendered Rough Surface Area 


19 
Therefore, at least M x M x N calculat 


As it can be seen from Figure 


TERCOM process. 


depends on the accuracy of the INS 


“M” 


Actually, it is obvious that 


the area considered should be larger ( 


quality is worse, 


thms (both MAD and MSD) are simp 


hand, TERCOM algori 


can be performed dur 


However, unnecessary calculations 


whole area is concerned. 


Next, the matrix cells are formed as seen in the figures above considering 


SSLM (short sample long matrix) map selection method for TERCOM [37]. The 


d sample size for the along track 


require 


simulations, it is selected nearly 10 km depending on the latitude of the area 
(Actually the numbers of the cells are taken to be constant.) Cross track errors 
depend on the accuracy of the INS. Depending on the typical 1 nmi/hr class INS 
which is generally used for cruise missiles, cross track sample size is selected 


approximately 1.85 km (1 nmi) assuming the worst case for cross track errors. 





Figure 20. Rendered Smooth Surface Area and Area with Unigueness 


On the other hand, velocity of the cruise missile is considered to be constant 
for the simulations moving from west to east direction (or vice versa) in order to 
investigate TERCOM concept. Before performing TERCOM algorithms, sigma-T 
and sigma-Z values are calculated in order to validate the roughness of the surfaces 


selected. 


First, considering perfect measurement and no error sources, selected 


profiles were determined. Here, the ambiguity term, é in equation (2.11) exactly 
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becomes zero as expected. Then, errors are added to the selected profile 


measurement values considering white noise. 


Monte Carlo simulations of 100 runs are performed for TERCOM 
simulations. Since, TERCOM is a batch process; true position fixes for navigation 
solutions are sought. Simulation results are given in Table 7 to Table 9 for different 


terrain types. 


Table 7. TERCOM Simulation Results for Rough Terrain 






































Initial INS Position Error (One axis, approximate) 400 m 
Height Measurement Standard Deviation (One sigma) 10 m 
TERCOM Map Grid Size (MxM) 21x21 
Number of Height Measurements for Correlation (N) 100 
Sigma-T of the Area Concerned 41.07 m 
Sigma-Z of the Area Concerned 1.88 m 
Correlation Method MAD MSD 
Percentage of False Fix 6% 9% 
Maximum False Fix Error (Total approximate error) 200 m 200 m 








Table 8. TERCOM Simulation Results for Smooth Terrain 























Initial INS Position Error (One axis, approximate) 400 m 
Height Measurement Standard Deviation (One sigma) 10m 
TERCOM Map Grid Size (MxM) 21x21 
Number of Height Measurements for Correlation (N) 100 
Sigma-T of the Area Concerned 722 m 
Sigma-Z of the Area Concerned 0.38 m 
Correlation Method MAD MSD 
Percentage of False Fix 100 % 100 % 
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Table 9. TERCOM Simulation Results for Terrain with Unigueness 


























Initial INS Position Error (One axis, approximate) 400 m 
Height Measurement Standard Deviation (One sigma) 10 m 
TERCOM Map Grid Size (MxM) 21x21 
Number of Height Measurements for Correlation (N) 100 
Sigma-T of the Area Concerned 46.4 
Sigma-Z of the Area Concerned 6.37 
Correlation Method MAD MSD 
Percentage of False Fix 1% 1% 
Maximum False Fix Error (Total approximate error) 150 m 150 m 


























From the results, it is seen that best position fix results are obtained with 
terrain with unigueness. However, it should be noted that the critical parameter is 
the sigma-Z value of the area concerned where standard deviation of the point-to- 
point changes in terrain elevation (i.e., the slope) are calculated. Eventhough, rough 
surface has larger sigma-T (standard deviation of height of the area) value, having a 
larger value of sigma-Z terrain with unigueness gives better correlation results than 
rough surface. Moreover, for smooth terrain, correlation algorithms do not give 


position fixes as expected. 


For the TERCOM process, several conclusions are achieved from the 


concept study and simulations performed. They are summarized as follows: 


1. Navigation solutions can be obtained for rough and unigue surfaces 


as expected. 


2. Correlation algorithm is simple but not smart. Many calculations 


should be performed in order to have a position fix. 
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It is thought that the algorithm was derived considering the 
capability of the computers of 1950's, performing only matrix 


calculations and simple mathematical operations. 


Physical meaning of MAD and MSD processes is the minimization 
of the area difference between the measured and the reference areas 


along the route of the missile. 


In the simulations, it was shown that MAD process shows better 
position fix than MSD process. For a terrain with small terrain height 
changes, MSD process neglects the small height difference terms and 
exaggerate the larger height difference terms. On the other hand, in 
MAD process absolute height difference terms are taken into account 


with same weights. 


The critical parameter for best terrain correlation is sigma-Z value of 
the area concerned where standard deviation of the point-to-point 
changes in terrain elevation (i.e. the slope) are calculated instead of 
sigma-T value where standard deviation of height of the area is 
calculated. In other words, the slopes of the area concerned are more 


critical than the roughness of the area for correlation. 


TERCOM process is independent of the target model where cruise 
missile is the target. Possible tracks for the missile are selected 
where tracks are the missile path formed by the terrain elevation file 
(DTED). Since, the target motion is not modeled, kinematical 


behavior of the system is not known. 


TERCOM process is actually a Maximum Likelihood Estimator 
(MLE) which uses “Least Sguares Estimation (LSE)” technigue. 
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Here, minimum error of the height measurements are sought for 


position fixes using least sguares (LS) estimation: 


(ah) tel) elk (2.16) 
£55 (k) =arg ex min) [20-40] [20-40] (2.17) 
where; 
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z(j): o Measurement of the terrain taken at “j 


h(j,x): DTED value of the related points with respect to taken 


measurement, z(/) 


X“(k): Minimum “MSD value times k” of the terrain height 


differences 


2.3. SITAN 
2.3.1. SITAN Fundamentals 


As it was stated in the first chapter, the major recursive TAN algorithm 
found in literature is SITAN which is proposed by Hostetler and Andreas [12]. In 
order to investigate SITAN in detail, first original work of Hostetler and Andreas 


[12] will be investigated in detail. 


The basic configuration for optimal terrain aided navigation is shown in 


Figure 21. This structure is typical of Kalman filtering in which nonlinear auxiliary 
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measurements are iteratively processed to estimate and compensate for the errors in 
a navigation system. At each measurement update time the current state estimate in 
conjunction with stored topographical data (i.e. terrain elevation data), is used to 
obtain a prediction of what the radar ground clearance measurement should be. The 
actual radar measurement is then compared with this predicted measurement, and 
their difference is processed by the Kalman filter to generate estimates of the 
navigation system's error states. The measurement matrix in this case is related to 
the downrange and cross range terrain slopes calculated from the stored data. The 
error estimates are then fed back to compensate the navigation system and thus 
provide an improved estimate of the actual state (position, velocity. etc.) of the 
system. This process is iterated many times e.g. every 30-50 m of distance traveled, 
as the system maneuvers along its trajectory, thus providing essentially continuous 


updating to the navigation system [12]. 





Suares 
Fir 
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Figure 21. SITAN Process [12] 


System equations for extended Kalman filtering (EKF) are derived 


considering navigation equations. True navigation state vector x is defined as [12]: 
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x = (2.18) 
Vy 
LY] 
where; 
x. Horizontal coordinates along eastward direction, 
y: Horizontal coordinates along northward direction, 
h: Height above sea level, 
v,: Velocity along x direction, 
v,: Velocity along y direction. 


Let x be the measured state vector for X from INS with the help of 


barometric altimeter, X the estimated state vector for X after updating, dx the 


optimal estimation of error vector 6x for x from the outputs of Kalman filter. 


For a constant sampling period 7, the recursion error state vector equation 


is as [12]: 
Ox(k +1) = O(k)- 6x(k) + W(k) (2.19) 


where; 
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ox =| Oh |: Error state vector, 








1 

0 
D(k)=|0 Transition matrix, 

0 

0 


SO O c r- o 
o OH ES 
SO — Sc o 3 
a Oo oy © 








Mk) =| WA) |: Process noise vector (White noises). 








In order to implement EKF, 1-D measurement h is needed, which is the 


difference between estimated relative height C,,, and measured relative height 


est ? 


Crease Cneas COMES from the measurement of radar altimeter; C,,, is the difference 


meas meas 


between estimated height above sea level from barometric altimeter (or INS), h, 


and terrain height horep from digital terrain elevation data based on the estimated 


position of (x,y) from INS. Thus öh is expressed by [12]: 
oh = Ca = C neas (2.20) 


where; 
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h Estimated height above sea level from barometric altimeter, 


baro * 


A A 


Aprep = he + Epren > 


A 


hprep: DTED (Terrain) height at the estimated position (x, Yy), 


h, Actual terrain height at the estimated position (X, f). 

meas 7 Padar * Wradar » 
haa: Radar altimeter measurement at the actual position (x, y), 
W „qa: Radar altimeter white noise measurement error. 


radar 


SITAN measurement process is shown in Figure 22. 


Now, expand terrain height difference measurement given in equation 


(2.20). 


öh = Cost > C = A (x, y) a (ho (x, y) + Erni] S [Prada (x, y) H Wradar ] (2:21) 


meas 
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Figure 22. SITAN Measurement Process 


Here, using Taylor series expansion, actual position is assumed to be near 


the estimated position. Therefore, 


uo (x, y) = aro (x, y) + ON jaro (2.22) 
Then, 
oh = m (x, y) T ON 9 | a. Mo (x, y) + Seen | ~ [A aaar (x, y) F Wadar | (2.23) 


The correlation between the estimated and the actual positions is the key 


point of the SITAN process [12]. Consider a fitted function f(x,y) to the terrain 


profile being expanded near (x, y) as shown in Figure 23. Then, 
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fan) J le Zy 5) 
X Oy 


where; 
of NEN 
h,= = : o Terrain slopes along eastward direction, 
x 
of , ye 
h, = a : Terrain slopes along northward direction, 
y 
Öx=X-x, 
Oy=y-y 


Here, terrain profiles h, and h, are needed for the EKF. 






hy 
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= , 
S fa.) = ' hz, y) 
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Horizontal Coordinates 


Figure 23. Terrain Stochastic Linearization (TSL) [12] 


TT 


(2.24) 


From the figure, following relations can be defined [12]: 


SOF)=f(G.y)+5f = f(xy) +h, x+h, öy (2.25) 
hy) = f (x,y) +m, (2.26) 
Io (X,¥) = f (X,y)+m, (2.27) 
where; 


m: TSL error at the actual point (x,y), 


m,: TSL error at the estimated point (x, y) 


Considering the actual position, following relation is valid [12]: 


[N (x, y) = hi (x, y) = hadar (x, y) = 0 (2.28) 


Then using eguations (2.23), (2.25), (2.26), and (2.27), 


oh=6h,,,,, —h, OEM, OV—(Eprep + W 


radar 


+m,+m,) 





(2.29) 


Wmeas 
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Here, terrain linearization errors m, and m, are included in the 
measurement error. Moreover, DTED error €,,,,, was modeled as white. Finally 


system measurement eguation in discrete form can be written as follows [12]: 
z(k) = H(k)-OX(k)+ Wea (K) (2.30) 
where; 


H (k) =|-h, x 1 0 o|: Measurement matrix, 


w,. (k): Measurement white noise 


Using linearization methods to the nonlinear system equations, the system is 
linearized. Therefore, EKF is implemented. Following part is the application of the 
standard EKF equations. EKF equations for SITAN are presented in Table 10 
considering standard EKF equations given in Gelb [63]. 


Implementation of the SITAN process is straight forward after the terrain 
slopes are modeled. However, the main problem of the process is the divergence of 
the KF. Due to highly nonlinear nature of terrain surfaces, filter divergence can 
occur especially when the linearization error is comparable to the measurement 
error. In these cases the standard EKF may yield unsatisfactory performance, and 
divergence can occur in which the actual estimation errors become orders of 
magnitude larger than the filter's own computation of their covariance [12]. Figure 
24 demonstrates this phenomenon for a simulation test case in which the initial 
position error standard deviations were 75 m and all other conditions were the same 


as in the prior simulation. 
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Table 10. Kalman Filter Eguations for SITAN Process 








SITAN System Error Model: 





OX(k AYDA) EZ) 
WA) N(0, O(k)) 


given in (2.19) 





System Noise Covariance Matrix: 





O(k) = Covi wk) wk) } 





(2.31) 





SITAN Measurement Model: 





2(k) = H(k)- OX (k) + w, (k) 
Wea (k) = NO, RK) 


given in (2.30) 





Measurement Noise Covariance Matrix: 





















































R(k) = COVE pegs (AY, KY} (2.32) 
Initial Conditions: 

OX, = N(6%,, P,) (2.33) 
Other Assumptions: 

E| WK) - Pos (K) | =O for all k (2.34) 
(Measurements are independent) 

State Estimate Propagation: 

Ox(k|k-1) = ®(k-1)-6x(k-1|k-1) (2.35) 
Error Covariance Propagation: 

P(k | k-1) = ®(k-1): P(k-1|k-1)-®(k-1)’ + O(k-1) (2.36) 
Gain Matrix: 

K(k) = P(k | k-1)-H(k)' -[H(k)- P(k| k-1)- H(k)" + R(k)] ' (2.37) 
State Estimate Update: 

Ox(k|k) = SX(k | k-1)+ K(k) -[z(k) - H(k)-Ox(k | k -1)] (2.38) 
Error Covariance Update: 

P(k|k)=[1—-K(k)-H(k)]-P(k | k-1) (2.39) 
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Figure 24. Standard EKF Divergence Problem [12] 


In order to improve EKF performance, modified stochastic linearization 
approach is used. However, single EKF for large errors actually can not perform 
good results. Therefore, parallel Kalman filters are used in order to estimate large 
position errors, especially large initial errors as shown in Figure 25. After initial 
errors are estimated within the accepted CEP values (i.e. ~30 m), single KF 


becomes sufficient for navigation purposes. 
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Figure 25. Parallel KF Configuration [12] 
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The selection of the convergent filter can be done guite easily by examining 
the residuals (i.e. estimated minus measured values of height) A, for each filter. A 
selection algorithm based upon the assumed whiteness property of the filter 


residuals that worked well in practice is to choose the filter with the smallest value 
of [12]; 





a: A, 
AWRS yy tiner = A 03 z (2.40) 


izl 


j'th filter 
where; 


AWRS n rec = Average Weighted Residual Squared of the j’th filter, 


J 


H,: Measurement vector containing the terrain slopes at the i'th time 
interval, 

P: Error covariance matrix, 

R,: | Measurement noise covariance matrix, 

N: Number of measurements processed, 


This AWRS value is the average weighted residual squared between the 
predicted ground clearance for each filter and the ground clearance measured by the 


radar altimeter for each time ¢,. The weighting factor is inherently calculated by 
each Kalman filter and is simply the expected variance of A, at each measurement. 


By examining the minimum AWRS values for each filter after a sufficiently large 
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number of measurements have been processed, the correct filter and its associated 


state error estimates can be chosen. 


2.3.2. Simulations and Discussion 


Simulations for SITAN are performed for both tracking and acguisition 
modes. In order to perform simulations, Simulink [58] is used. Mathematical 
models described in the previous section are used for trajectory and INS models in 
order to obtain 1.0 nm/hr INS guality by adding white noise terms to horizontal 
positions, altitude and horizontal velocities. Terrain slopes are derived considering 


the gradients of the height values of the related DTED files. 


For the simulations, three special terrain types are selected: 


1. Rough terrain, 


2. Smooth terrain, 


3. Mountainous terrain. 


Some properties of these selected terrains for TAN are given in Table 11. It 
should be noted that, these properties satisfy terrain reguirements for the 


simulations. 


Simulation model details will be presented in the following section. In this 
section, SITAN characteristic simulation results will be presented. First, horizontal 
position errors for tracking mode are performed for three different terrain types. 


Simulation parameters for tracking mode are given in Table 12. Here, it should be 
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noted that initial position error is less than the grid size of the DTED considered 


(i.e. less than 100 meters for DTED Level 1). 


Table 11. Terrain Parameters for SITAN Simulations 




















Terrain Type Rough Smooth Mountainous 
Mean height of the terrain profile 1093 m 1104 m 1177 m 
Sigma-T 71.9 m 34.1 m 212.9m 
Sigma-Z 16.3 m 3.7m 23.1m 

X, 670.2 m 1309 m 1302 m 





























Table 12. SITAN Simulation Parameters for Tracking Mode 


























Initial INS position deviation (one axis) 80 m 
Initial vehicle velocity 240 m/s 
Initial INS east velocity bias 0.5 m/s 
Initial INS north velocity bias 0.5 m/s 

INS horizontal position standard deviation 5m 

INS altitude position standard deviation 3m 

Radar altimeter standard deviation 3m 

INS velocity standard deviation 0.3 m/s 




















SITAN filter works at 1 Hz. In other words, it gives updates at every 1 
second. Simulations are performed for 100 seconds of operation time. In actual 
systems, INS is updated recursively considering SITAN position corrections. 
Hence, INS errors become zero at discrete SITAN updates. However, in the 


simulations, in order to show SITAN characteristics, INS error model is not 
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updated; and only first 100 seconds of operation is considered. Simulation results of 


tracking mode for different terrain types are shown from Figure 26 to Figure 31. 


Northward Position Error vs. Time 
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Figure 26. Rough Terrain Northward Position Error vs. Time 
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Figure 27. Rough Terrain Eastward Position Error vs. Time 
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Northward Position Error vs. Time 
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Figure 28. Smooth Terrain Northward Position Error vs. Time 
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Figure 29. Smooth Terrain Fastward Position Error vs. Time 
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Northward Position Error vs. Time 
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Figure 30. Mountainous Terrain Northward Position Error vs. Time 
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Figure 31. Mountainous Terrain Eastward Position Error vs. Time 
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As it can be seen from the figures above, SITAN improves position errors 
for rough and mountainous terrain types. On the other hand, due to slope 
determination process in SITAN, solutions have serious jumps for mountainous 
terrain type. This can be explained by the severe slope changes in the mountainous 
terrain modeling. As a result of this, SITAN works better for rough terrains. 
However, by TSL or other linearization methods as explained in the original paper 
[12], navigation solutions can be improved. In fact, these technigues depend on the 


terrain selected; and, extra work is reguired for terrain linearization. 


Next, SITAN simulations for acquisition mode are performed. Here, initial 
position error is assumed to be greater than DTED grid size. 25 parallel KF’s are 


66 3 99 
1 


used in the simulations as shown in Table 13. Here, the index indicates the 
related grid for initial position. For example, if the initial position error was 2” short 
along longitude and 2” long along latitude considering INS outputs, actual initial 


position would be at i=5 where i=13 was the INS index. 


Table 13. Parallel KF Structure for SITAN Acguisition Mode 









































AZ hays +2*3" i=5 i=10 i=15 i= 20 i= 25 
A= Ans +3" i=4 i=9 i=14 i=19 i=24 

A= Ans i=3 i=8 ERİ i=18 i=23 
A= Ams 73" i=2 i=7 i=12 i=17 i=22 
A= Ay 2*3" i=l i=6 i=11 i=16 i=21 




















Time:t=k M= Ms —2*3" U= Ms 3" M= ys HF Hys +3" M= Hg +2*3" 


Note: Index i=13 gives h(k) at position (Å ys » Hys ) Of INS, at time “t = t”. 


Index i=1 gives dh(k) at position (4,,—2*3", H,,, —2*3") of INS for DTED Level 


1, at time k. 
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Hence, the simulations are performed for acguisition mode. Simulations are 
performed only for rough terrain in order to show acguisition performance. 
Simulation parameters for acguisition mode are given in Table 14. Here, initial 
position errors are given for both axes in order to determine initial position index. 


From the simulations, determination of the initial position index is reguired. 


Table 14. SITAN Simulation Parameters for Acguisition Mode 





























Initial INS position deviation (northward axis) -200 m 
Initial INS position deviation (northward axis) -180 m 
Initial position index (according to Table 13) 25 

Initial vehicle velocity 240 m/s 

Initial INS east velocity bias 0.5 m/s 

Initial INS north velocity bias 0.5 m/s 
INS horizontal position standard deviation 5m 
INS altitude position standard deviation 3m 
Radar altimeter standard deviation 3m 

INS velocity standard deviation 0.3 m/s 




















Simulations are performed for rough terrain for both minimum AWRS filter 
and the central filter. Horizontal position errors are given in Figure 32 and Figure 
33. Minimum AWRS filter index versus time is given in Figure 34. Here, central 
SITAN filter results are also presented in order to show filter divergence. From the 
simulations, it can be seen that in order to obtain correct navigation solutions for 
large initial position errors, parallel KF’s should be used. Moreover, from Figure 34 
initial position index obtained is exactly the same with the simulation initial 
condition which means that the correct initial position is found from the parallel 


filter structure SITAN simulations. 
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Northward Position Error vs. Time 


e deği, 


2 
a 
Ww 
= 
x 
O: 
2 
O 
© 


[aau] 


INS 
— — Central SITAN KF 


Min. AWRS SITAN KF 





[second] 


Figure 32. Northward Position Error vs. Time for Acquisition Mode 
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Figure 33. Eastward Position Error vs. Time for Acguisition Mode 
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KF İndex vs. Time 
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Figure 34. Minimum AWRS KF Index vs. Time for Acguisition Mode 


For the SITAN process, several conclusions are achieved from the concept 


study and simulations performed. They are summarized as follows: 


1. 


SITAN is a recursive TAN technigue which uses EKF unlike 
TERCOM which is a batch process. 


SITAN performance depends on the linearization of the terrain 
profiles since terrain slopes are reguired for the KF measurements. 
For large position errors, divergence can occur due to linearization 
errors in the EKF. In order to get rid of this, modified terrain 


linearization technigues and parallel KF structure are used. 


SITAN improves position errors for rough and mountainous terrain 
types. However, due to slope determination process in SITAN, 
solutions have sometimes serious jumps for mountainous terrain 


type. This can be explained by the severe slope changes in the 
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mountainous terrain modeling. Therefore, linearization of the terrain 
profiles is very critical especially for mountainous terrains in 


SITAN. 


SITAN performance is better than both INS and terrain grids unlike 
TERCOM. In TERCOM, error can not be better than the terrain grid 


dimensions. 


SITAN performs better for smaller position errors due to terrain 
linearization. Due to this fact, for large initial position errors 


TERCOM or SITAN with parallel KF structure must be used. 


SITAN is a tracking process (i.e. it tracks the actual path with 
minimum errors) where TERCOM is an acguisition process (i.e. it 


estimates the initial position of the target). 


VATAN 


2.4.1. VATAN Fundamentals 


As it was stated in the first chapter, one of the interesting TAN algorithms 


found in literature is VATAN which is proposed by Enss and Morrell [42]. In order 


to investigate VATAN in detail, first original work of Enss and Morrell [42] will be 


investigated in detail. 


VATAN is a recursive TAN technigue which uses Viterbi Algorithm (VA). 


VA is a maximum a posteriori (MAP) estimator that estimates a seguence of system 


states from a seguence of observation values [42]. Viterbi algorithm is actually a 


dynamic programming technigue for estimation which uses past information. 
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The state and observation seguences are denoted by [42]: 


X = (Xos X, ) (2.41) 
Z = (Busey Zi) (2.42) 
where; 


x(k 
X, -| ( | Vehicle’s position at time /,, 


Z, = [2(k)] : Measured terrain elevation at time /,. 


The VA consists of the computation of a metric function L, that is a 
measure of the likelihood of each state value being the true state at time k; L, can 
be computed recursively using conditional probability density functions p(x, | X,) 
and p(z,|x,) as follows [42] based on the assumption that the system dynamics 


are Markov; that is, the state at k +1 is conditionally independent, given the state at 


time k, of the state at any previous time: 


k k 
L,(%,) = max p In p(X, |.) +11 p(X,)+ 2 In p(z, | 9) (2.43) 
Ipari izl i=l 
k 
= In p(z, |X,)+ max È In p(x, |.) + hh] (2.44) 
kili 
where; 
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L,(x,) =In p(x,): Initial condition of Z,. 


The optimal estimate X, is that X, for which Z, is maximum. The value 
X, , that maximizes equation (2.44) for each x, is termed the survivor. Denoted 
S,(Xx,), it is used to generate MAP state sequence estimates. The MAP state 
sequence estimate xX=(x,,...,%,) can be generated via the following recursive 


procedure [42]: 


X, =argmax(L,(7,)) oa 
Kd =S, x.) 
Ea = 1104) (5) 


% 2) = GOY EE) (2.46) 


The recursion in equation (2.44) is a filter, providing state estimates based 
on the system dynamics and observations. For an observable linear system model 
with Gaussian noises, equation (2.44) is functionally equivalent to a Kalman filter 
and equation (2.46) is equivalent to a fixed interval smoother (e.g., a Rauch-Tung- 
Striebel smoother). These equivalences suggest that the VA is a suitable 


replacement in applications that use Kalman filtering [42]. 


For the TAN problem, the VA has two significant advantages over the EKF 
used in SITAN [42]: 
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1. Metric function is computed for all possible state values makes the 
VA much more robust than the EKF in situations where the 
observations do not strongly support a single estimate of the state 


value. 


2. The nonlinear relationship between vehicle position and measured 
terrain elevation can be represented exactly with the VA but must be 


approximated for the EKF. 


In VATAN, the VA generates optimal MAP vehicle position estimates using 
the terrain elevation beneath the vehicle as its observation. In order to implement 
VATAN, conditional observation and; state transition densities in eguation (2.44) 


are needed as well as an initial value of the metric L,.In the original VATAN 


paper [42], very simple models are used to obtain the reguired densities. 
Parameters reguired for the VATAN technigue are given as follows [42]: 


Nominal terrain height (Actual terrain height with zero measurement errors): 


Pre ONI Nays OI — hadar % ) (2.47) 


Measured terrain elevation (Observation used in VATAN): 


Z, = SRK) = Ig, (Xp) = hs (Xp) — Prater He) (2.48) 
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The measurements are assumed to be unbiased, independent and Gaussian. 


Therefore; E[hns(X,)] = hws) and Eh ml huu (%,) with variances 


radar 


o; (k) and o; (k).Then [42]: 


Elz, | x] = Es Eid GIN = Mer E) (2.49) 


o (k) = 05 (to; (6) (2.50) 


Thus the conditional observation probability density function is: 





(2.51) 


Male — e) 
ee Proti) 207 (6) 


The state transition density p(x,,,|x,) describes the states’ evolution with 
time. Given a known velocity vector x, that is constant over the 7 second sample 


interval from /, to ¢,,,, 


the state’s evolution is [42]: 


Xi İT (2.52) 


Since the vehicle’s velocity is provided by the INS, it is not known 


precisely. This uncertainty is dealt with by modeling the INS velocity as a Gaussian 


INS 


random variable with mean x and variance Cx , indicative of the INS 
k 
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precision. Thus the random variable x,,, conditioned on Xx, is a random variable 


with mean and variance [42]: 


Elx, ]=% + ED, “TT (2.53) 


a = Ons T? (2.54) 


Since x, is assumed to be Gaussian, by equation (2.52), x, is Gaussian for 


k>0. Moreover, the metric function L, is initialized considering equation (2.43) 


by [42]: 


L) — In p(x) = In p(x,"*) (2.55) 


Here, it should be noted that VATAN models are derived considering a 
simplified INS model. Actually, INS error model used in SITAN can also be used 
for VATAN implementation. 


2.4.2. Simulations and Discussion 


Simulation results for VATAN are presented from the original work of Enss 
and Morrell [42]. They performed simulations for VATAN using four different 


terrain types: 


1. Typical rough terrain, 
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2. Flatterrain, 


3. Mountainous terrain, 


4. Sloped and flat terrain. 


Flight paths regarding the terrain types are shown in Figure 35. 





Northing (100s of meters) 
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Figure 35. Contour Plot for VATAN Simulation Terrain Types [42] 


Simulations are performed for horizontal position errors (mean errors 
performed with Monte Carlo simulations and deviation errors simulated in tracking 
mode) and the results are compared with SITAN. Simulation results for different 
terrain types show that VATAN consistently performs as well as or better than 


SITAN implementation. VATAN performs as well as SITAN in moderately rough, 
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sloped terrain and it exceeds SITAN’s performance in very flat or very rough 


terrains [42]. 


However, there exist some drawbacks of the VATAN technigue. VATAN's 
major limitation is the increased computational capacity necessary to implement the 
VA when compared to an EKF [42]. In the original paper, the two-dimensional VA 
has only been implemented in a discrete state space; a continuous state-space 
implementation of the two-dimensional VA would improve VATAN's accuracy 
and could result in a substantial reduction of the computational capacity necessary 


to implement VATAN [42]. 


For the VATAN process, several conclusions are achieved from the 
investigation of the original paper of Enss and Morrell [42]. They are summarized 


as follows: 


1. VATAN is a recursive TAN technigue like SITAN. However, since 
the past measurements are stored and used its performance is said to 


be better than SITAN. 


2. From the paper, it is shown that VATAN performs better results for 


all terrain conditions (both very rough and flat terrains). 


3. VATAN uses VA which is a maximum a posteriori (MAP) estimator 
that estimates a seguence of system states from a seguence of 
observation values. VA is actually a dynamic programming 


technigue for estimation which uses past information. 


4. The major disadvantage of VATAN is the limitation of the increased 
computational capacity necessary to implement the process. 


Actually, VA is also used as a radar tracking algorithm. By 


99 


investigating this paper, implementation of modern radar tracking 


algorithms to TAN has been inspired. 


In this chapter, major TAN algorithms have been investigated in detail with 
their fundamentals described in original references and the simulations performed. 
Simulation model details are not presented in this chapter; since, they will be given 
in the following chapter. Several conclusions have been obtained from the detailed 


study of the major algorithms and they have been discussed in the chapter. 
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CHAPTER 3 


IMPLEMENTATION OF TARGET TRACKING ALGORITHMS 
TO TERRAIN AIDED NAVIGATION 


In this chapter, implementation of target tracking algorithms to TAN is 
presented. First, general information about modern target tracking algorithms are 
given. Next, PDAF and TSF data association algorithms and their general 
implementations are investigated. Then, PDAF and TSF implementations to TAN 
are presented. At the end of the chapter, a simple simulation model is developed for 
the mid-course flight of the cruise missile. Finally, simulations are performed with 
the implemented TAN algorithms and the results are compared with the major TAN 


methods. 


3.1. Target Tracking Background 


In the first section of the chapter, a historical background about target 
tracking will be introduced. Major developments in multi-target tracking over the 
past four decades and how algorithms developed primarily for tracking air targets 
will be discussed. Then, target state estimation algorithms like Kalman filtering and 


association algorithms fundamentals will be investigated. 
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Eventhough tracking problems can be found in many applications, e.g., 
ocean surveillance and submarine tracking, most tracking algorithms have been 


developed for air targets [64]. 


A tracking problem is defined by the targets of interest, the sensors that 
collect the measurements, and the environment in which the targets move and 
sensors observe the targets. The basic functions in multi-target tracking consist of 
prediction, association, and estimation and they are shown in Figure 36. When 
measurements are received, the current tracks are predicted to the time of the 
measurements and associated with the measurements. Then the associated 
measurements are used to update the state estimates of the tracks. Although these 
functions are not always performed seguentially, they are present in most tracking 


algorithms [64]. 





Figure 36. Basic Tracking Functions [64] 


Prediction and estimation are single target state estimation functions in the 
absence of measurement uncertainty. Prediction difficulty depends on target 
dynamics and sensor revisit time. On the other hand, when the origins of the 
measurements are uncertain, e.g., when clutter or multiple targets are present, the 
measurements have to be associated with other measurements or tracks before the 


target state estimates can be generated. Therefore, data association establishes 
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tracking as a separate discipline from traditional state estimation in target tracking 


[64]. 


Chong, et al [64] investigated target estimation and data association 
algorithms in detail for ground target tracking. However, they discussed the subject 
from the historical point of view and investigated the algorithms with the related 
references in detail. In the following sections, target estimation and data association 
algorithms will be summarized considering the helpful reference of Chong, et al 


[64]. 


3.1.1. Target State Estimation 


Target state estimation is an important component of any multi-target 
tracking algorithm. The association of measurements to tracks reguires the 
prediction of the target state of each track to the time of the measurements so that 
the measurement to track likelihood can be computed. Accurate state prediction is a 
key to good association performance. Once an association decision has been made, 
the output of the tracker consists of updated state estimates of the tracks using the 


associated measurements [64]. 


Target state algorithms can be grouped according to the algorithms applied 


as follows: 


1. Linear Estimation Algorithms: 


These algorithms assume linear target motion and observation 
models and provide estimates of the target state by means of linear 


transformations [64]. 
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2. 


a. Alpha-Beta-Gamma Filters: 


These constant coefficient filters estimate the target position and 
velocity from position measurements only. The alpha-beta filter 
assumes a second order model driven by white noise for the target 
dynamics while the alpha-beta-gamma filter assumes a third order 
model. In either case, the filters can be considered as steady state 


Kalman filter [64]. 


b. Kalman Filter: 


The Kalman filter has been the standard approach to filtering for 
linear systems since its development in the earlier sixties [64]. 
Details of Kalman filtering have been discussed in several chapters 


of the study. 


Adaptive Filters: 


When a target maneuvers, the model no longer matches the 


dynamics and performance will degrade. Several approaches have been 


developed to detect maneuvers and adapt the filter to the target dynamics 


in real-time [64]. 


a. Parameter Adjustment: 


The structure of the filter is fixed. However, the filter will 
monitor its own performance (such as the size of the residuals) and 
adapt parameters (such as the process noise covariance or the 


Kalman filter gain) when a target maneuver is detected [64]. 
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3. 


b. State Augmentation: 


This approach uses different dynamics models when a maneuver 
has been detected. For example, before maneuver, a constant 
velocity model is used. When a maneuver has been detected, the 
filter switches to an acceleration model with higher state dimension 
and switches back to the original model when the maneuver is 


determined to have ended [64]. 


Multiple Models: 


When the measurement does not contain sufficient information, an 


incorrect decision may be made, resulting in poor performance. 


Therefore, algorithms that maintain multiple target dynamic models 


have been developed. These algorithms compute the probability of each 


model being true given the measurements and generate a target state 


estimate as a weighted sum of the estimates given the individual models 


[64]. 


a. Static Multiple Models: 


These models assume that the true target motion model is static 
and contained in a fixed set of models. Because the target model 
does not change with time, this approach is not appropriate for 


maneuvering targets [64]. 


b. Model Sequence Pruning: 


The optimal multiple model estimator requires a filter for each 
possible model sequence hypothesis. Since the number of model 


sequences and thus the number of filters increases exponentially with 
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time, the optimal estimator is not practical. An obvious sub-optimal 
approach is to prune the least likely model seguences according to 


their probabilities (64). 


c. Generalized Pseudo Bayesian Estimator: 


The Generalized Pseudo Bayesian (GPB) method is a suboptimal 
approach that reduces the number of filters by merging model 
seguences that end up with the same fixed length sub-seguences 


[64]. 


d. Interacting Multiple Models: 


The Interacting Multiple Model (IMM) algorithm is one of the 
most popular algorithms for tracking maneuvering targets because of 
its relatively simple implementation and its ability to handle 


complicated dynamics [64]. 


e. Variable Structure Interacting Multiple Models: 


While IMM has been successfully used in several applications, 
having a fixed model set has its disadvantages. Variable Structure 
Interacting Multiple Model (VSIMM) approach is used to track 
ground targets moving over roads and open field. The target motion 


models reflect the mobility of a target for different conditions [64]. 


4. Nonlinear Estimation: 


Many dynamic models or observation models do not satisfy the 
linear assumptions. Therefore, approaches for estimating the state of 


nonlinear systems have been developed. 
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a. Extended Kalman Filter (EKF): 


The non-linearity of the dynamic and observation models can be 
linearized about a nominal trajectory, and then a Kalman filter can 


be developed with the linearized model which is called EKF [64]. 


b. Gaussian Sum Approximations: 


The EKF assumes that the conditional probability distribution 
can be approximated reasonably accurately by a Gaussian 
distribution. When this approximation is not valid, the conditional 
probability distribution of the states given the cumulative 
measurements can be approximated by a sum of Gaussian 


distribution [64]. 


c. Nonlinear Filtering: 


This optimal nonlinear filtering algorithm has nice features such 
as the ability to update the probability distribution of the states due 
to non-detections. However, implementation is computationally 
intensive since it reguires discretization of the state space and 
performing the integration by a summation. Thus, even though the 
algorithm has been known for many years, it has seldom been used 


[64]. 
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3.1.2. Data Association 


When the origins of the measurements are uncertain, e.g., when clutter or 
multiple targets are present, the measurements have to be associated with other 
measurements or tracks before the target state estimates can be generated. 
Association is what distinguishes target tracking from traditional state estimation 


and establishes tracking as a separate discipline [64]. 


Data association algorithms can be classified according to whether they 
focus on single targets or consider explicitly the presence of multiple targets and 
whether association decisions are made using single or multiple scans of data. The 
early algorithms tend to focus on single scan and single targets, while the recent 
algorithms deal with multiple scans of data and multiple targets. In general, 
algorithms that consider multiple targets and use multiple scans of data perform 


better but reguire more computations [64]. 


Data association algorithms can be grouped according to the algorithms 


applied as follows: 


1. Single Target Track Formation: 


These track formation algorithms initiate tracks from seguences of 


measurements without considering competition from other tracks [64]. 


a. “M” out of “N” Test: 


A track is tentatively initiated from a single measurement. A 
validation gate is then established around this measurement and a 
measurement falling inside this gate becomes part of the track. When 


there are “M” detections out of “N” scans of measurements, then the 
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track is formed or confirmed. This method is very simple but does 


not provide a score on the confidence of the track [64]. 


b. Likelihood (Ratio) Test: 


In the likelihood tests, tracks are declared as confirmed (or 
deleted) when the likelihood or ratio exceeds (or falls below) a 


certain threshold [64]. 


2. Single Target Track Maintenance: 


These algorithms associate measurements with the existing tracks 
without considering the presence of other tracks. Thus a measurement 


may be associated with multiple tracks [64]. 


a. Nearest Neighbor: 


In this method, the measurement that is closest (according to 
some distance measure) to the track is associated with the track from 
the multiple measurements. This approach makes a hard decision 
based on a single scan and is very easy to implement. However, it 


does not perform well in high density situations [64]. 


b. Track Splitting: 


This is basically applying the likelihood function (or ratio) 
approach to track maintenance. For every measurement that falls in 
the validation gate, the track is split. Each track is scored using a 
likelihood function as discussed before. The track is pruned when the 
likelihood falls below a threshold. This approach makes soft 


decisions based upon multiple scans of data. Because of its 
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3. 


computational requirements and limited performance, this approach 


is no longer popular [64]. 


c. Probabilistic Data Association (PDA): 


Instead of associating a single measurement with a track, this 
approach probabilistically associates all measurements in the 
validation gate. The PDAF is an all-neighbors association algorithm. 
It is fairly easy to implement and has been shown to perform better 


than the nearest neighbor approach in high clutter [64]. 


d. Optimal Bayesian Approach: 


The PDAF is a suboptimal approach since the association event 
only considers the current measurements. On the other hand, the 
optimal Bayesian approach will consider all possible association 


hypotheses up to the current time [64]. 


Multiple Target Track Maintenance: 


Association performance can be improved when the algorithms 


consider explicitly the presence of multiple targets and recognize that a 


single measurement cannot belong to multiple tracks [64]. 


a. Optimal Assignment: 


The optimal assignment approach, also sometimes called global 
nearest neighbor, is the coordinated version of nearest neighbor. 
Instead of selecting the measurement that is closest to a track, this 
approach selects the set of measurements that is closest to the set of 


tracks according to some global distance measure subject to the 
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4. 


constraint that two tracks do not share a single measurement, and 


two measurements do not appear in the same track [64]. 


b. Joint Probabilistic Data Association (JPDA): 


This is the extension of PDA to multiple targets. The tracks for a 
known number of targets are assumed to have been initiated and the 


problem is to associate the measurements to the tracks [64]. 


Multiple Scan Coordinated Association: 


Both the measurement and the target motion models have 


uncertainty. Therefore, the single scan decisions may not be the correct 


associations. Thus association performance can be improved by using 


multiple scans of data. The core of all multiple scan algorithms is the 


evaluation of track likelihoods, which can be used for both track 


formation and maintenance. Thus multiple scan algorithms generally can 


be used for both track formation and association [64]. 


a. Integer Programming: 


This approach was the fist multiple scan algorithm and integer 
programming problem can be solved by branch and bound or other 
methods. However, this algorithm was improved to other multiple 


scan algorithms such as multiple hypothesis tracking [64]. 


b. Multiple Hypothesis Tracking (MHT): 


Multiple hypothesis tracking delays making hard decisions when 
there is not sufficient information to make a good decision. 


Alternative hypotheses are formed to represent the ambiguities and 
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5. 


each hypothesis is evaluated. The MHT is conceptually simple but 
computationally intensive since the number of hypotheses grows 


exponentially [64]. 


c. Multi-Dimensional Assignment: 


Traditional MHT requires the explicit expansion and evaluation 
of many hypotheses. Successful implementation requires the use of 
sophisticated hypothesis management techniques to handle the 
combinations. During the last decade, alternative optimization based 
methods that do not require the explicit expansion and evaluation of 
hypotheses have been developed. Such algorithms are easier to 


implement and computationally more efficient [64]. 


Tracking Without Data Association: 


Several approaches have been proposed to perform tracking without 


an explicit association function. Instead of dealing with individual target 


states and individual measurements, these approaches treat all targets 


and measurements as components of one system, and estimate the 


system state directly without explicitly forming association hypotheses 


[64]. 


a. Symmetric Measurement Equations: 


In this approach the original measurements on the targets are 
converted into a new set of measurements that are symmetric 


functions of the original measurements [64]. 
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b. Multi-target Nonlinear Filtering: 


The individual target motion and measurement models can be 
aggregated into a multi-target motion model given by the conditional 
probabilities and a measurement model given by the likelihood 
function. Then, the same nonlinear filtering method developed for a 
single target can be used (conceptually at least) for tracking multiple 


targets [64]. 


As it can be seen from the historical point of view, target tracking is a 
comprehensive subject. In this study, implementation of some of these algorithms to 
TAN is done. As a result of this, some of the algorithms summarized above will be 
discussed in detail in the following sections. Then, they will be implemented for 


TAN applications. 


3.2. Probabilistic Data Association Filter (PDAF) 


3.2.1. Theory 


The PDA algorithm calculates in real-time the probability that each 
validated measurement is attributable to the target of interest. This probabilistic 
(Bayesian) information is used in a tracking filter, the PDA filter (PDAF) which 


accounts for the measurement origin uncertainty [52]. 


The following assumptions are made to obtain the recursive PDAF state 


estimator (tracker) [52]: 
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1. 


There is only one target of interest whose state evolves according to 


a dynamic eguation driven by process noise. 
The track has been initialized. 


The past information about the target is summarized approximately 


by; 

p| x(k)] Z“ | =N[ x(k); 8(k | kK), P(k|k-1)]z (3.1) 
where, 

N [x(k); X(k | k -1), P(k | k - 1)] : Normal probability density function. 
x(k): Argument, 

x(k|k-1): Mean, 

P(k|k-1): Covariance matrix. 


At each time, a validation region is set up. 


. Among the possibly several validated measurements, at most one of 


them can be target-originated, if the target was detected and the 


corresponding measurement fell into the validation region. 


The remaining measurements are assumed to be false alarms or 
clutter and are modeled as independent identically distributed 


measurements with uniform spatial distribution. 
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7. The target detections occur independently over time with known 


probability. 


These assumptions enable a state estimation scheme to be obtained, which is 


almost as simple as the Kalman filter, but much more effective in clutter [52]. 


The probabilistic data association algorithm associates all valid observations 


with a track. For each validated observation, an updated estimate x,(k|k) is 
computed. A probability of correct association 4, is computed for each such track. 


Then a combined track is formed from the weighted average of these tracks: 


KAİ) >) B,-X,(k | k) (3.2) 


For multiple targets, the same process occurs although the probability 
calculations are more complex which is called Joint Probabilistic Data Association 


Filter (JPDAF) [65]. 


In Figure 37, general PDAF implementation is shown. Then, the following 


approach can be implemented in order to perform PDAF algorithm [65]: 
1. The set of validated measurements is computed. 
2. For each validated measurement an updated track is computed. 


3. For each updated track an association probability Ø, is computed. 


The calculation of this probability can be quite complex and 


dependent on the assumed clutter densities. However, it is normally 
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adequate to set Ø, proportional to the normalized innovation for the 


association. 


A 
x(klk) 





Figure 37. PDAF Implementation [65] 


4. A combined (average) track is computed. 


5. A combined average covariance can also be computed although this 


can become guite complex. 
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6. A single prediction of the combined track is then made to the next 


scan time and the process is repeated. 


The PDAF and JPDAF methods are appropriate in situations where there is a 
high degree of clutter. The great advantage with the PDAF method is that you are 


never wrong. The problem is you are also never right [65]. 


PDA procedure is summarized in Figure 38 and detailed PDAF eguations 


are given in Appendix. 


State estimate at tk-1 State prediction State covariance at tk-1 
x(k - 1)k -1) x(kIk - 1) P(k -4|k -1) 


N 


Predicted measurement State prediction covariance 
Z(klk-1) P (kik -1) 


Innovation calculation and Innovation covariance 
measurement validation Sik) 
i= mik) 









Measurements 


Filter gain 
wik) 










Evaluate association 
probabilities — %.(k) 


v(k) 


Updated state estimate 
x(klk) 















Effect of measurement 
origin uncertainty — Pik) 









Updated state covariance 
P(klk) 






Figure 38. PDAF Procedure [66] 
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3.2.2. Implementation of PDAF to TAN 


3.2.2.1. Implementation Methods of PDAF for TAN 


Terrain Aided Navigation (TAN) algorithms estimate the position of a 
moving vehicle by comparing the measured terrain profile under the vehicle to a 
stored elevation map. Therefore, the critical part of the TAN is the terrain elevation 
database, i.e. DTED for the common military applications. Batch and recursive 
algorithms are used for TAN as explained before. For batch algorithm, i.e. 
TERCOM, only DTED is used in order to estimate the position of the vehicle. On 
the other hand, for recursive algorithms, like SITAN, VATAN, etc. besides DTFD, 


the system dynamics should also be modeled. 


Application of the target tracking algorithms to navigation problems have 
been investigated in several papers given in Oingtang, et al [40], Dezert [43] and 
Maksarov and Durrant-Whyte [67]. In the paper of Oingtang, et al [40], TAN using 
PDAF was investigated for the batch algorithm. Association probabilities have been 
derived using the MSD function of the TERCOM process and performance of the 
TAN using PDA and TERCOM has been compared. In the paper of Dezert [43], 
PDAF has been used in order to improve the accuracy of a strapdown INS using 
landmark detections. Maksarov and Durrant-Whyte [67] used multiple hypothesis 
technigue (MHT) algorithm for an autonomous mobile vehicle with multiple sonar 


sensors for range measurements. 


As it can be seen from the papers investigated, Oingtang, et al [40] is 
directly related with the Ph.D. study. The main difference of the Ph.D. study from 
Oingtang, et al [40] is the real-time application of PDAF to TAN which will be 


discussed in detail. 


First application of a target tracking algorithm to TAN problem should be 


discussed. Consider the assumptions to obtain a recursive PDAF estimator again: 
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1. There is only one target of interest whose state evolves according to 


a dynamic eguation driven by process noise. 


e Target of interest is the cruise missile and its motion can be 
modeled. In MHT and JPDA multiple targets are 


considered with more complex algorithms. 


2. The track has been initialized. 


e Initial conditions of the vehicle’s motion can be modeled. 


3. The past information about the target is summarized. 


e Cruise missile dynamics model gives information about its 


motion. 


4. At each time, a validation region is set up. 


e CEP of the vehicle is determined by the quality of the INS 
used. The validation region for the measurements is set 
considering 36 position error bound of the INS horizontal 


positions. 


5. Among the possibly several validated measurements, at most one of 
them can be target-originated, if the target was detected and the 


corresponding measurement fell into the validation region. 


e DTED area considering the 36 position error bound of the 
vehicle can be used considering only one of the height 


measurements rely on the target. 
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6. The remaining measurements are assumed to be false alarms or 
clutter and are modeled as independent identically distributed 


measurements with uniform spatial distribution. 


e False measurements are modeled considering the height 
channel covariances due to radar and INS height 


measurements. 


7. The target detections occur independently over time with known 


probability densities. 


e Radar height measurements can be modeled as time 


independent occurrences. 


From the assumptions of PDAF, it is seen that PDAF can be applied to 
TAN. Data association problem comes from the DTED height differences used for 
TAN. It is known that one of the grid of the DTED batch considered gives the 
correct position of the vehicle. Now, the implementation of PDAF algorithm to 


TAN can be discussed. 


First, consider other TAN algorithms used for both batch and recursive 


algorithms. In Figure 39, batch algorithm concept is summarized. 


In the batch algorithm, only the height measurements and their relations with 
the related DTED are considered. For a period of time, measurements are taken and 
correlation can be obtained using TERCOM process which is actually a maximum 
likelihood estimator. For the batch algorithm, the model of the vehicle motion is not 
reguired. Therefore, procedure is simple and larger DTED area should be used for 
the calculations. As a result of this, it can be concluded that batch algorithm can be 
successfully used in cruise missiles with considerably accurate INS. In fact, INS 


guality used in cruise missiles is around 1 nm/hr. Then, a few position updates 
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during the mid-course phase of the operation with the help of TERCOM algorithm 


will be sufficient for the navigation solution. 


Corrected INS State 
< DTED Grid Size 





DTED Grid 
(Larger Size) Measured Height 
Difference (0/,) = 
t=, 
Actual Path SZ a 15 
Batch Solution 
INS Path 


t=t INS Height Difference (hys ) 


Figure 39. Batch Algorithm for TAN Solution (Acquisition Mode) 


In recursive algorithms, TAN solution is done continuously. In order to 
achieve a complete navigation solution, the motion of the vehicle should also be 


modeled. In Figure 40, recursive algorithm for TAN concept is shown. 
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Corrected INS State 
<< DTED Grid Size 


Measured Height 
Difference (0h, ) 
DTED Grid 
(Smaller Size) 
| J . et | t=k+2 
Actual Path A J 


Recursive Solution 
INS Path t=k at each time step 


INS Height Difference (hys ) 


Figure 40. Recursive Algorithm for TAN Solution (Tracking Mode) 


As it was discussed in earlier chapters, there are various recursive TAN 
algorithms proposed in literature. The well-known recursive algorithm is SITAN. 
Recursive TAN algorithms generally use estimation theory in order to solve the 
navigation problem. In SITAN, extended Kalman filtering, where in VATAN, 
Viterbi algorithm (a maximum a posteriori state sequence estimator) is used. 
Various recursive algorithms are also proposed as using maximum a posteriori 


estimation and optimal Bayesian estimation. 


Due to the nonlinear dynamics of the navigation system using terrain 


information, algorithms for recursive TAN solutions require generally complex 
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calculations. In SITAN, linearization is done using EKF for the terrain. Several 


disadvantages of the recursive TAN algorithms can be summarized as follows: 


1. TAN reguires terrain information for the navigation solution and the 


dynamics of the system is highly nonlinear. 


2. Eguations derived for recursive algorithms are generally complex 


and needs considerable calculation work. 


3. Real-time application for the TAN solution is generally impractical 
for high velocity vehicles like cruise missile due nonlinear 


characteristics of the system. 


4. In SITAN, terrain linearization and terrain slopes are reguired in 
order to apply extended Kalman filter eguations which are actually 


critical stages for TAN solution. 


In this study, a new recursive TAN algorithm is investigated, which uses 
PDA, a target tracking algorithm. In Figure 41, real-time PDAF application for 
TAN is shown. Consider the PDA approach given in the previous section for TAN 


again: 


1. The set of validated measurements is computed. 


e Measurement gate is taken as 30 position error bound of 
the INS and the invalid possibilities for the height 
differences Oh (k) are discarded. 


2. For each validated measurement an updated track is computed. 
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e Here, INS error model is used due to its linear 
characteristics. The updated tracks Ox, (k) are computed 


for all valid points in the gate using PDA filtering. 


3. For each updated track an association probability 4, is computed. 


e Using height differences for each valid element of the grid, 


association probabilities 5; are calculated. 


4. A combined (average) track is computed. 


5. A combined average covariance can also be computed although this 


can become quite complex. 


6. A single prediction of the combined track is then made to the next 


scan time and the process is repeated. 


e Estimated error state dx(k|k) is computed considering 


equation (3.2). Then the estimated error state becomes: 


dx(k|k) = XP, EZİN) (3.3) 


Next, derivation of the PDAF equations for TAN is done. PDAF equations 
for TAN are implemented considering standard real-time PDAF equations which 


are given in Appendix. 
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Actual Height 
Difference (0/,) 






5% (k | k) 


t=k+1 


INS Height Difference (hys ) 


Figure 41. PDAF Implementation for TAN Solution 


3.2.2.2. PDAF Eguations Implemented for TAN 
3.2.2.2.1. Past Measurement Information 


In the TAN algorithm, the only measurement is the height difference Oh 
given in eguation (2.23) which was given for SITAN in the previous chapter. 
However, since the 30 horizontal error bound of the vehicle is estimated 
considering the guality of the INS used, a batch of height differences around the 
INS position can be obtained. Rewrite SITAN eguations considering the barometric 
height given as the INS height. Then, barometric height definition can be rewritten 


as follows: 
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hys (x, Y) = Haro Y) (3.4) 


Since DTED are used in the simulations where the heights are given with 
respect to related longitudes and latitudes, instead of using eastward and northward 


positions (x, y), longitudes and latitudes (4,4) are selected in the equations. Then, 


estimated and trajectory (real) positions can be defined as: 


(3, J) = (His > Ays) (3.5) 
(X,Y) = (Maj Aira) (3.6) 
where; 

(x,y): Estimated eastward and northward position, 


(Mms> Ais): Longitude and latitude of the INS position, 


(x,y): Trajectory eastward and northward position, 


(raj Ayaş): Longitude and latitude of the trajectory position. 


Considering the measurements to be taken at discrete time steps k, SITAN 


measurement eguations derived according to Figure 22 can be rewritten as follows: 
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hys (k) = hys (Las (k), Ans (AD) 


= horeo (Hins (k), Ans (K) + wus (6) + Cn (R) Pn 
hrs (E) = horeo Hin O) Ary OD + Wir (K) + Cren (K) (3.8) 
A aaar (K) = C meas (K) + agar (K) (3.9) 
Oh(K) = C peas K) — Cost (k) (3.10) 
SICK) = [Apren (ins (K), Ans OD) + Wins İİ N 


z [ Peis (Maj (4), Apa (K))+ W adar (4) | 


Now, consider the DTED batch for each height difference which is shown in 
Figure 42. Then, for each grid node, eguation (3.11) can be written as follows 


considering the grid index i: 


öh (k) = [pre Uk) 4;(k)) + Wrys(k)] 
-| horen (Hra (K), Aaj OD + adar CK) | 


(3.12) 
i=1,...,m(k) 
where; 


i Index of the DTED grid node, 


m(k): DTED grid size (selected as square of an odd number for INS 


position to be at the center of the DTED grid) 
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Here, it should be noted that INS height difference is at i=13 fora 5x5 
DTED grid size. 


30 s ~ 500m 


INS Position 


(5x5 DTEDI Grid) 
A=A,,-2"3" 





H= Hpg —2*3" H= His H= Hys +2*3" 
Figure 42. Height Differences Batch Used in PDAF 


Height difference at the INS position can be found for a DTED grid size of 
m(k) as: 


Öh nys (k) = Oh, (k) (3.13) 


zim) 


From Figure 42, it can be seen that the only actual measurement is 


h, = hys at time step k. Other 6h,, i—l,....m(k) values are derived using the 
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defined DTED grid. Using the DTED heights given at the positions around the 
given INS position, height difference batch at time step k is obtained. 


For the TAN algorithm, past height difference measurements are averaged in 
order to smooth the effects of the past measurements. In order to achieve this, 
height measurements are put in a buffer and then, the average is used as the 
measurement batch to the PDA filter. For the study, buffer size of 20 to 30 (i.e., 20- 
30 seconds of measurements) was sufficient. As the new measurements come, the 


oldest measurements are eliminated. 


Consideration of the past measurement information can be summarized as 


follows: 


s+k 


öh, (s+4)="-Y"6h(k) s>k>0 (3.14) 
ie s 4 


where; 
öh, : Average of the height difference at position related to index i, 


S: Buffer size. 


In the same manner, height difference batch matrix can be written as: 








[ oh, (s+k) 
oh, (s+k) 
Sh (stk) =| Oh, (sth) < N. s>k20 (3.15) 
Omya, (SR) 
[pe oe e b Yİ ii 
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3.2.2.2.2. Measurement Validation 


In order to determine validation region for measurements, consider standard 
PDAF eguations in the Appendix. In the TAN algorithm, validation region given in 
equation (A.3) is directly taken as the 30 error bound of the vehicle. However, 
height difference measurements in the batch matrix must be valid for the 
calculations. In the eguation considered, measurements are the height differences 
oh,,i =1,...,m(k) . Innovation covariance S(k) given in equation (A.4) contains the 
system height state covariance P(k|k—-1), and the radar measurement noise 


covariance matrix R(k). Therefore: 


P(k | k =) = Oi ;R(k) = Cy (3. 16) 


Then, equation (A.3) becomes: 


Sh (k) [Pk k- +R EN) <y 
oh, (ON = y (oi. + Orr) 


h(k) < 


rlo tou) (3.17) 








Gate threshold y is taken as 16 (40 error bound) considering 99.9989% of 


the measurements to be in the gate as in the original reference of Kirubarajan and 


BarShalom [52]. Height difference measurement values, óh (k), which are not 
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valid according to eguation (3.17) are eliminated by assigning very large values of 


öh; (k). Hence, in the data association process for the measurements given in 
eguation (A.33), probability of the invalid measurements become zero (i.e. e =0) 


and have no effect in the TAN solution. 


3.2.2.2.3. State & Covariance Estimation, Update and Prediction 


For the TAN algorithm, generic PDAF eguations are directly used using a 


modification for the definition of the states. Filter eguations are given as follows: 


INS Error Model: 
öx(k +1) = O(k)- dx(k)+ W(k) (3.18) 
where; 


Ox(k) = [Or ;6rE;öh; öyN; vE IN : Navigation error states vector, 


örN: Northward position error state, 
örE: Eastward position error state, 
oh: Height position error state, 
OvN: Northward velocity error state, 


OvE: Eastward velocity error state, 
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D(k)= State transition matrix with sample time T. 


SOS Som 
S-OO [S 
S: © es SS 
S = o 53 
— O CON o 








WK) = [15,3 (K), Won (K), Wop (K), Woy (K), Ws (k)|: INS error state white 
noises where w(k)=N(0,07) with mean 


zero and variance o; to the related state. 


PDAF Measurement Model: 


z (k) = H p (E) OX) + W peas E) (3.19) 


where; 


z;(k)=öh, (k): Average of the height difference at position related to 


index i, 


H (k) = [0 0 1 0 0]: Height measurement matrix, 


m 


w,, (k)=N(0,0? 


meas radar 


): Measurement White Noise with mean zero and 


variance o° 


radar * 


After defining system and measurement models, PDAF eguations can be 


implemented for TAN. PDAF eguations contain Kalman filter eguations with 
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association probabilities. For PDAF eguations, innovation form of Kalman filter 


eguations are used. 


Propagation (Prediction) Eguations: 


Ox(k|k-1) = ®(k-1)-6x(k-1|k-1) (3.20) 
P(k|k-1) =@(k-1): P(k-1| k-1)-®(k-1)" + O(k-1) (3.21) 
where; 


P(k): State covariance matrix, 








o 0 0 0 0 
0 a, 0 0 0 
O(k) =| 0 sä 0 0 |: System noise covariance matrix 
0 0 0 œ 0 
0 0 0 0 07] 


State and Covariance Update: 


S(k) = H, (k): P(k|k-1)-H,(k)" + R(k) (3.22) 
K(k) = P(k|k-1)-H,(k)" -S(k)" (3.23) 
öX(k | k) = OX(k|k-1)+K(k)-v, (k) (3.24) 
where; 
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S(k): Innovation covariance matrix, 


K(k): PDA filter gain, 


1 0 0 0 0 
H,(k)=|0 1 0 0 0 PDAF measurement matrix, 
001 0 0 
m(k) 
vp(k) = > BA) ve (k): Combined PDAF innovation, 
i=0 
vp (k) = z, (k) 2, (k|k-D): PDAF innovation states, 
Zp(k)= | orn, pOrk,; oh, ] s PDAF measurement states, 


Z2p(k|k-1)=H,(k)- öx(k|k-1): Measurement state estimation, 


Gn 0 0 
R(k)=| 0 o; 0 |: PDAF measurement noise covariance matrix 
0 0 o; 


radar 


Here, it should be noted that PDAF measurement states z, (k) are different 
from the actual measurement states z,(k). Since, only height difference 
measurements are taken into account, position updates are not available with the 
height measurement matrix H, (k). As a result of this, PDAF measurement matrix 


H,(k) is defined in order to make position corrections with the same filter gains of 


height corrections. 
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PDAF measurement states z,(k) are defined such that northward and 


eastward positions are calculated according to the index i as shown in Figure 42. 
Moreover, due to position inaccuracies along horizontal coordinates, horizontal 


position white noises o, and o. terms are added in the measurement noise 
covariance matrix. Height difference measurements z,(k) are used for the 
determination of the conditional probability of the event 8.(&) instead of z,(k); 
since, positions are not actually measured. With the use of z,(k), filter gains 


obtained for height difference states are directly used for position states. 


For the position error definitions, consider Figure 42 again. As it was given 
in equation (3.13), INS position index is at ins = ln) +] for a DTED grid size 
of m(k) where the grid is selected as a square. At the INS position index is, 


OrN ws =0 and rE ws <0; since no position correction exists for the INS position. 


Using geometrical relations for the DTED grid index and horizontal positions, 


following definitions can be done for horizontal position errors: 


ses) ia cil : | Ota (3.25) 





Jmk) 2 





ob, = e á ron), (3.26) 


where; 


örN,;: o Northward position error at index i, 


örE,: Eastward position error at index i, 
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m(k): DTED grid size, 


ceil(X): Function which rounds the element of X to the nearest integer 


towards infinity. 


d,: DTED spacing along latitude direction, 


y 


d.: DTED spacing along longitude direction. 


x 


Horizontal position error definitions are shown in Figure 43. 





INS Position 





H= Lys —2*3" H= Ens H= Lys +2*3" 


Figure 43. Horizontal Position Error Definitions 
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As an example, consider the horizontal position errors at index i=2 for 5x5 
DTED grid size as shown in Figure 43. Using eguations (3.25) and (3.26), 


horizontal position errors can be found as follows: 


ön, =[ 2445. el iz) D aaa, 








örE, , ea za ee 


As it can be seen from the calculated results, horizontal position errors, 
which are determined from equations (3.25) and (3.26), can be directly used for INS 
position updates. Therefore, by determining the index of the position from TAN 


algorithms, INS error model can be updated. Moreover, as it was stated earlier, 
horizontal position white noises oc^, and o. terms are added in the measurement 


noise covariance matrix in order to model position inaccuracies along horizontal 


coordinates. 


For the state covariance update, equations (A.15) to (A.17) are used. 
However, the conditional probability of the false events is zero (i.e. 8,(k4) = 0) for 
the TAN application which will be explained in the following section. Moreover, 
for the spread of the innovations term P(k), only the measured state z,(k) will be 


considered. Hence: 


P(k|k)=P(k|k-1)-K(k)-S(k)-K(k)" + Pk) (3.27) 
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m(k) 


Pk) = «|S B(k)-v(k)-v (ky! 9-10" | Ky (3.28) 


where; 


m(k) 


v(k)= X B(k)-v,(k): Combined height difference innovation, 
i=0 


v,(k) = z,(k)—2,(k | k—1) 


v;(k)=öh, (k)-H(k)- x(k |k —1) : Height difference innovation states 


Here, it should be noted that P(k) term will be effected only from the height 


channel. Therefore, in state covariance matrix P(k|k), P(k) term will be added 


only to the height channel. 


3.2.2.2.4. The Probabilistic Data Association 


Association probabilities are calculated considering the height differences 


used in parametric PDA eguation (A.32). In this eguation, determination of P, and 
P, parameters are critical. Probability of detection of a target originated 
measurement P, must be one; since the height difference measurements grid is 


formed virtually from the related DTED within the 30 horizontal error bound of the 
vehicle position. If no measurements are taken, then the height measurement grid 


could not be formed. Probability of measurements in the gate P, is also taken one 
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considering the used DTED grid is the 36 horizontal error bound of the vehicle 


obtained from the guality of the INS. Actually measurement gating process 


eliminates the impossible height difference solutions according to eguation (3.17) 


derived. Therefore: 


P, =1,P, =1 


D 








Then from equation (A.34), b 








e 





0 and eguation (A.32) becomes: 


Pk) = n i=1,...,m(k) 
ey 
j=l 
where; 
-Lv (T-S v; (k) feb. 
eze? : Given in equation (A.33) 


$ 


3.2.2.2.5. Summary of PDAF Equations for TAN 


(3.29) 


(3.30) 


PDAF equations derived in the previous sections for real-time TAN 


application is summarized in Table 15. Initialization of Kalman filters is generally 


done by setting state covariance matrix as a coefficient of system covariance matrix: 


P, =a -Olk), 


(typically æ =10) [65] 
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(3.31) 


Table 15. PDAF Eguations for TAN Process 








INS Error Model: 





Kk AYDA) EZ) 
mk) = N(0,0(k)) 


given in(3.18) 





System Noise Covariance Matrix: 





O(k) = Coviw(k)w(k)"} 





given in (3.21) 





PDAF Measurement Model: 





z; (K) = H, (k) ÖZE) + css (k) 


m meas 


Wreas (k) = N(0, O padar) 


given in (3.19) 





Association Probabilities: 





€; 


BAK) =) aw i=l,...,m(k) 


İ 





JA 





given in (3.30) 





Initial Conditions: 











5X, = N(6%,,P)) (3.32) 
Other Assumptions: 
E| W(k)- Wg. (k) | =O for all k (3.33) 


(Measurements are independent) 





State Estimate Propagation: 





Ox(k|k-1) = ®(k-1)-6x(k-1|k-1) 


given in (3.20) 





Error Covariance Propagation: 





P(k| k-1)=®(k-1)- P(k-1|k-1)-®(k-1)! + O(k-1) 


given in (3.21) 





PDAF Gain Matrix: 





K(k) = P(k|k-1)-H,(k)" -S 


given in (3.23) 





State Estimate Update: 











Ox(k | k) = 8x(k | k-1)+ K(k)-v,(k) 





given in (3.24) 
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Table 15. PDAF Fguations for TAN Process (Continued) 








Error Covariance Update: 


P(k|k)=P(k|k-1)—-K(k)-S(k)-K(k)" + P(k) given in (3.27) 








m(k) 


E - ; i given in (3.28) 
P(k)= Kof BA) v;(k)-v;(k) —v(k)-v(k) | <) 























3.2.2.3. Discussion of Real-time PDAF Implementation for TAN 


In this section, implemented PDA method for TAN will be discussed. The 


advantages of the PDA approach for TAN solution can be summarized as follows: 


1. Real-time TAN solution can be obtained with a single PDA filter. 


2. PDA filter can be used for both batch and recursive TAN solution. 
For batch solution, larger grid size is selected for navigation solution. 
For recursive solution, horizontal positions are calculated recursively 


in relatively small DTED grids. 


3. Since past measurements are taken into account, smoothing of the 


measurements in the filter is achieved which decreases errors. 


4. Since INS error model is used for navigation solution, application of 


the filter is simple and the filter is linear. 


5. Batch size of the DTED area concerned can be changed. Both larger 
DTED areas for acquisition mode or smaller DTED areas for 


tracking modes can be selected using the same filter. 
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6. 


Results of the filter are good for both recursive and batch algorithms. 
The results compared with SITAN and TERCOM algorithms will be 


discussed in the simulations section. 


The difference of the PDA approach from Oingtang, et al [40] is also 


summarized as follows: 


1. 


In the paper of Oingtang, et al [40], TAN using PDAF was 
investigated for the batch algorithm. The motion of the vehicle is not 


modeled. 


In the paper of Oingtang, et al [40], the batch algorithm obtained 
using PDAF actually uses maximum likelihood approach as used for 


TERCOM. Therefore, association probabilities Ø, are calculated 


with the help of the MSD function used in TERCOM. 


In the paper of Qingtang, et al [40], performance of the TAN using 
PDA and TERCOM has been compared. It is stated that PDA was 
used in order to improve the performance of TAN compared to 


TERCOM. 


In the Ph.D. study, real-time PDAF implementation is done. By 
using the error model of the INS used in the vehicle, system 
dynamics is modeled. Using PDAF, error states of the system are 


estimated. 


In the Ph.D. study, PDAF equations are directly implemented for the 
TAN solution. Association probabilities obtained from height 
difference measurements for each element of the DTED grid 
concerned are used for position updates considering the index of the 


DTED grid. 
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Simulation results of the implemented PDA filter will be presented at the 


“Simulations” section of the chapter. 


3.3. Multiple Hypothesis Tracking (MHT) and Track Splitting Filter 
(TSF) 


3.3.1. Theory 


In classical multiple-target tracking, the problem is divided into two steps, 
association and estimation. Step 1 associates contacts with targets. Step 2 uses the 
contacts associated with each target to produce an estimate of that target's state. 
Complications arise when there is more than one reasonable way to associate 
contacts with targets. The classical approach to this problem is to form association 
hypotheses and to use MHT. In this approach, alternative hypotheses are formed to 
explain the source of the observations. Each hypothesis assigns observations to 
targets or false alarms. For each hypothesis, MHT computes the probability that it is 
correct. This is also the probability that the target state estimates that result from 
this hypothesis are correct. Most MHT algorithms display only the estimates of 
target state associated with the highest probability hypothesis [68]. 


The model used for the MHT problem is a generalization of the recursion for 
general multiple-hypothesis tracking. This recursion applies to problems that are 
nonlinear and non-Gaussian as well as to standard linear Gaussian situations. In this 
general case, the distributions on target state may fail to be independent of one 
another (even when conditioned on an association hypothesis) and may reguire a 
joint state space representation. This recursion includes a conceptually simple 


Bayesian method of computing association probabilities [68]. 


143 


Numerous books and articles on multiple-target tracking examine in detail 
the many variations and approaches to MHT problem. Many of these discuss the 
practical aspects of implementing multiple target trackers and compare approaches. 
In addition to the full or classical MHT as defined by “Reid” and “Mori et al.”, a 
number of approximations are in common use for finding solutions to tracking 
problems. Examples include joint probabilistic data association and probabilistic 


MHT [68]. 


Multiple hypotheses tracking (MHT) is a deferred decision logic in which 
alternative data association hypotheses are formed whenever there are observation 
to track conflict situations. Then, rather than combining these hypotheses, as in the 
JPDA method, the hypotheses are propagated in anticipation that subseguent data 


will resolve the uncertainty (69). 


The original MHT method, denoted Reid’s algorithm, was first presented by 
Reid [70]. There are two basic approaches to MHT implementation. The first 
(hypothesis-oriented) approach follows the original work of Reid [70]. It maintains 
the hypothesis structure from scan to scan and continually expands and cuts back 
(prunes) the hypotheses as new data are received. At each scan, a set of hypotheses 
will be carried over from the previous scan and composed of one or more tracks that 
are compatible with all other tracks in the hypothesis. Compatible tracks are defined 
to be tracks that do not share any common observations. Then, on the receipt of new 
data, each hypothesis is expanded into a set of new hypotheses by considering all 
observation-to-track assignments for the tracks within the hypothesis. Again, as new 
hypotheses are formed, the compatibility constraint for tracks within a hypothesis is 


maintained [69]. 


An alternative (track-oriented) approach [71] does not maintain hypotheses 
from scan to scan. The tracks formed on each scan are reformed into hypotheses 
and the tracks that survive pruning are predicted to the next scan where the process 


continues [69]. 
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In Figure 44 and Figure 45, the operations of MHT that are reguired by both 


implementation methods are summarized. 


Observations 







From New 
Tracks / Hypotheses 






Surviving 
Tracks / Hypotheses 











Track 
Prediction 


Hypothesis / Track 
Evaluation / Deletion 










Most Likely 
Hypotheses 


User 


Figure 44. MHT Logic Overview [69] 


New 
observation set — «— s" f 
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correlate with this observation i 
> 


| 
Start now Merge involved 
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Finished? == o | 


Yes No 





Loop clusters o ————— 


New cluster? 
No Yes 


Generate feasible 
corretation hypotheses 


Evaluate probability 
of each hypothesis 


Combine sımilar tracks and hypotheses 
Prune unlikely hypotheses 


Finished? ooo. 


Yes | No 


Filter all tracks 
Delete bad tracks 


Split clusters if possible 
Dele 


te Hii M cistes 


Figure 45. High-level Flow Chart of MHT Algorithm [69] 
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The Multiple Hypothesis Tracking (MHT) filter maintains separate tracks 
for each possible associated observation. At each time step, the predicted 
observation is used to establish a validation gate and for each measurement that is 
found in this validation gate, a new hypothesis track is generated. Thus a single 
track is split into “n” tracks, one associated with each valid measurement, plus one 


track (usually denoted 0) for the no-association hypothesis [65]. 


Fach of these new tracks is then treated independently and used to generate 
new predictions for the next time step. Since the number of branches into which the 
track is split can grow exponentially, the likelihood function of each split track is 
computed and the unlikely ones are discarded. The MHT algorithm works on 


complete seguences of observations [65]. 


In MHTF, every validated observation z,(k) is used to establish a new 
track, x,(k|k). In addition the “false alarm” and/or “missed observation” 
hypothesis also generates a track, X,(£ |k). These tracks are propagated forward to 


the next gate and again each track is associated with each valid observation, 


z,(k+1) and the tracks are again split into tracks associated with each possible 
pair-wise association, X, (£+1|£+1). Probabilities or likelihoods, 4,, of correct 


track histories are maintained to prune the resulting hypothesis tree [65]. 


In Figure 46, MHTF implementation is shown. 
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Figure 46. MHTF Implementation [65] 


Then, the following approach can be implemented in order to perform 
MHTF algorithm [65]: 


1. A predicted observation and validation gate are computed. 


2. All validated observations are associated plus no-association 


hypothesis. 


3. The track is updated separately with each validated hypothesis. 
4. A likelihood associated with correct association is computed. 
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5. The likelihood of the each entire track sequence is computed. 


6. Some pruning of the hypothesis tree may take place. 


7. Each track hypothesis is now independently predicted forward to the 


next time-step. 


8. Producing as many new tracks as associated measurements (plus no 


track solution). 


9. The process repeats. 


There are three points to note about the MHT and TSF algorithm [65]: 


1. A unity detection probability (no missing data) is assumed. 


2. Likelihood pruning method does not work well with long 
measurement sequences as it becomes dominated by old 
measurements. One “hack” around this is to use a fading-memory 


window. 


3. Method is dominated by computational and memory requirements of 


the splitting algorithm. 


MHT algorithm is good in situation with low clutter rates but high track 
uncertainty (crossing tracks, maneuvering targets, etc). Practically, the algorithm is 


dominated by the approach used for pruning unlikely target hypotheses [65]. 
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TSF is proposed by Smith and Buechler [72] and older than the original 
MHT method presented by Reid [70]. In TSF, a tree of hypotheses is kept for each 
target individually, and a maximum likelihood criterion is used to prune the tree. On 
the other hand, Reid's MHT constructs a tree of all possible hypotheses, including 
all possible new track initiations at every time step. Reid discusses a number of 
strategies to prune the tree in order to achieve reasonable computation times. In the 
Ph.D. study, TSF is implemented for TAN due to INS error model characteristics. 
Since, horizontal INS error bound is estimated for the cruise missile and errors do 
not change rapidly, implementation of TSF for TAN became sufficient for 


navigation solution. 


3.3.2. Implementation of TSF to TAN 


As it was stated in the previous section, TSF is an older method than the 
original MHT method. TSF is a recursive branching algorithm for multiple-object 
discrimination and tracking consists of a bank of parallel filters of the Kalman form, 
each of which estimates a trajectory associated with a certain selected measurement 
seguence. The measurement seguences processed by the algorithm are restricted to 
a tractable number by combining similar trajectory estimates, by excluding unlikely 
measurement/ state associations, and by deleting unlikely trajectory estimates. The 
measurement seguence selection is accomplished by threshold tests based on the 


innovations seguence and state estimates of each filter [72]. 


TSF and MHT methods are similar except for hypotheses formation. 
Consider the TSF and MHT approach given in the previous section for TAN again 
[65]: 


1. A predicted observation and validation gate are computed. 
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e Measurement gate is taken as the 30 horizontal error bound 
of the INS and the invalid possibilities for the height 
differences Oh (k) are discarded. 


2. All validated observations are associated plus no-association 


hypothesis. 


e Every grid position (i.e. index) in the 36 horizontal error 
bound of the INS is considered to be one of the possible 


navigation solutions. 


3. The track is updated separately with each validated hypothesis. 


4. A likelihood associated with correct association is computed. 


5. The likelihood of the each entire track sequence is computed. 


e Navigation solution is assumed to be one of the grid 
positions in the 30 horizontal error bound of the INS. 
According to the index of the grid position, there exist “n x 
n” possible tracks (i.e. hypothesis) for each time step where 
“n x n” denotes the batch size of the DTED considered. The 


likelihood of the each possible track sequence is computed. 


6. Some pruning of the hypothesis tree may take place. 


e Number of possible tracks is limited considering INS error 
characteristics. According the small position error changes 
of the INS for small periods of time where the TAN 
algorithm is applied, it is assumed that possible tracks are 


in the 36 horizontal error bound of the INS where each 
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track follows the grid position from INS position to all grid 
positions in the INS horizontal error bound. Then, tracks 
with minimum likelihoods are selected for the navigation 


solution. Hence, hypotheses are pruned. 


7. Each track hypothesis is now independently predicted forward to the 


next time-step. 


e Navigation solution is found for each possible track using 
standard Kalman filter equations as given in the reference 
papers for MHT/ TSF procedure. Using a definite number 
of minimum likelihood values of the entire track sequences, 


navigation solution is achieved. 


Actual Height 
Difference ( Ôh, ) 






e. 7 â (klk) 
K 


Compare “4 ” for each track. 
(For all NxN tracks) 
Search for definite number of 


minimum “A” values. 
INS Height Difference (hys ) 


Figure 47. MHTF Implementation for TAN for a Single Time Step 
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Details of the TSF procedure applied for TAN and derivation of the TSF 
eguations are discussed in detail in the following section. In Figure 47, TSF 


application for TAN for a single time step is shown. 


3.3.3. TSF Eguations Implemented for TAN 


TSF maintains separate tracks for each possible associated observation. At 
each time step, the predicted observation is used to establish a validation gate and 
for each measurement that is found in this validation gate, a new hypothesis track is 
generated. Thus a single track is split in to “n” tracks, one associated with each 
valid measurement, plus one track (usually denoted 0) for the no-association 


hypothesis [65]. 


Fach of these new tracks is then treated independently and used to generate 
new predictions for the next time step. Since the number of branches into which the 
track is split can grow exponentially, the likelihood function of each split track is 


computed and the unlikely ones are discarded [65]. 
The TSF procedure works as follows [65]: 
1. TheTSF algorithm works on complete seguences of observations. 


2. The probability that a given branch seguence of observations (from 


root to leaf) is correct. 


3. Thel'th seguence of measurements up to time k: 


Z" = (a, Dz, ()) (3.34) 
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©” is the event that the sequence Z” is a correct track. 


. Then the likelihood function for this event is clearly: 


A(O) = P(Z* | O") = Plz, (0,2, (189) (3.35) 


Z* the cumulative set of all measurements up to time k: 


k 
MO") = (le, |Z") (3.36) 
z 


. Linear and Gaussian distribution for the likelihood function is 


assumed: 
k,l l E Tru a cer i 
MO" )=c, -exp = 2 GS (7) v) (3.37) 
j=l 
where, 


v(j)=z(j)-Z(j| 7-1): Innovation between track and measurement 


. Modified log-likelihood function is defined as: 


A(0*”) 


k 





Ak) = 20g) İst) (3.38) 


. Modified log-likelihood function is recursively computed from: 


Mk) = A(k-l) + v" (k) S (kK) vk) (3.39) 
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10. Each track is updated using standard Kalman filter eguations. 


11. A “goodness of fit” and test for accepting a track is that A(k)<d. 


12. Definite numbers of tracks are accepted for each time step. 


As it can be seen from the TSF procedure, for each track, standard Kalman 
filter eguations are used. Hence, PDAF state estimation, state and covariance update 
and prediction eguations which are given in the previous sections are used in order 


to include horizontal position inaccuracies to the navigation solution. 


The critical part of the TSF procedure implemented for TAN is the track 


formation and track pruning steps. These steps are summarized as follows: 


1. Navigation solution is assumed to be one of the grid index followed 
by some of the tracks in the 36 horizontal error bound of the INS. 
According to the index of the grid position in the 36 horizontal error 
bound, there exist “n x n” possible navigation solutions where “n x 


n” denotes the batch size of the DTED considered. 


2. Possible tracks are different from grid indices. At the initial time 
step, there exist “n x n” possible tracks from INS position grid to all 
possible grid positions as shown in Figure 48. The modified log- 
likelihood of the each possible track seguence is computed from 


eguation (3.38) as follows: 
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AO") 


k 


A(k) = 2) |- 2, DSTO) v 0) 


for i=1..."nxn 
where, 


m(k): DTED Grid Size (DTED grid is taken as a square.) 


Corrected INS State 


Measured Height 
Difference ( h, ) 


; kf 
DTED Grid Ap Td K 






ye 
t=k+3 
A, t=k+2 
1. Compute all tracks from INS 
position to all grids. 
2. Select best “M” tracks. 


3. Compute all tracks from best 
INS Path i=k “M” tracks. 


4. Repeat procedure. 
INS Height Difference (hys ) 


Actual Path 


Figure 48. TSF Track Formation and Pruning 


3. At the following time step, best “M” tracks of the existing “n x n” 
possible tracks are selected. Hence, hypothesis is pruned. In order to 


apply pruning, several methods can be applied. Using 4,(k)<d for 
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accepting tracks or selecting minimum 4(k) values for i—l...M 


can be applied. Then, “n x n” possible tracks from “M” accepted 
tracks to all possible grid positions are selected. Modified log- 
likelihood function is recursively computed from eguation (3.39) for 


each formed new track as follows: 


Ayk =A,(k)+ vy, (AS (kv, (k+1) 
fori=1..M, j=1..."nxn" 


. Procedure defined at step 3 is done recursively in order to obtain 


navigation solution. 


. For real-time navigation solution, accepted tracks can be used in 
several methods. Selecting and using the results of the best track 
which gives minimum likelihood is actually a dynamic programming 
method which was discussed for Viterbi algorithm. On the other 
hand, using mean value of the selected tracks? results can also be 
used which is actually a kind of data association process discussed 
for PDA algorithm. However, from the simulations which will be 
discussed in the following sections, all tracks converge to the same 
index of the DTED grid. Actually, for sufficiently rough surfaces this 
is the expected result of the TSF. 


. In order to decrease the effects of the old measurements, modified 


log-likelihood function defined in eguation (3.39) can be used by a 


weighting factor as follows: 


A, (k+1) = Kyp Ak) i (£+1)-S'(k+1)- vy, (k +1) Aa) 
forizl..M, j=1..."nxn" l 
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TSF eguations derived for TAN application is summarized in Table 16. 
Initialization of Kalman filters is again done by setting state covariance matrix asa 


coefficient of system covariance matrix given in eguation (3.31). 


As it can be seen from Table 16, standard Kalman filter eguations are used 
for each track considering horizontal position error inaccuracies. Track formation, 


pruning and real-time solutions are also summarized. 


Table 16. TSF Eguations for TAN Process 








INS Error Model: 
öx(k +1) = O(k)- dx(k)+ Wk) given in(3.18) 


w(k) = N(0, OK) 


System Noise Covariance Matrix: 














O(k) = Coviw(k)w(k)"} given in (3.21) 
TSF Measurement Model: 
Z.(k)= H, (k)-0x,(k)+Ww,,.,.(k) given in (3.19) 


w, (k)=N(0,07 


meas radar ) 





Measurement Noise Covariance Matrix: 




















R(K) = Cov{W peas (K)W peas (K) 3 given in (3.30) 
Initial Conditions: 

OX, = N(6%,, P) given in (3.32) 
Other Assumptions: 

E| WK): mass (k) | =O for all k given in (3.33) 


(Measurements are independent) 





State Estimate Propagation (for each track): 














Ox(k|k-1) = ®(k—1)-6x(k-1|k-1) given in (3.20) 
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Table 16. TSF Eguations for TAN Process (Continued) 








Error Covariance Propagation (for each track): 





P(k|k-1)=&(k-1)-P(k-1]/k-1)-0(k-1)" + 0(k-1) 


given in (3.21) 





TSF Gain Matrix (for each track): 





K(k) = P(k|k-1)-H,(k)' -S(k)' 


given in (3.23) 





State Estimate Update (for each track): 





Ox(k | k) =8x(k | k-1) + K(k)-v,(k) 


given in (3.24) 





Error Covariance Update (for each track): 
P(k|k)=[1-K(k)-H(k)]-P(k|k-1) 


given in (2.39) 





Modified Log-likelihood Function (Initial Track Formation): 








(k) = zog MO k 


Dv (DSV, O) 


Ck 


for i=1..."nxn" 


given in (3.38) 





Recursively Computed Modified Log-likelihood Function: 





A (k+1) = Kyr A (k) + Vir, (k +1): STk +D vy, (k +1) 
fori=1..M, j=1.."nxn" 





given in (3.39) 





Track Pruning: 




















Select best “M” tracks from “M x n x n” tracks such that; (3.41) 
A(k+T)=min A,(k*T) = forij=1..M 
State Estimate Propagation (for all tracks): 
dx(k |k-l)= average| 5x,(k |k- D| for i=1...M 3.42) 








3.3.4. Discussion of TSF Implementation for TAN 


In this section, implemented TSF method for TAN will be discussed. The 


advantages of the TSF approach for TAN solution can be summarized as follows: 


1. Real-time TAN solution can be obtained with parallel TSF’s. 
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. Application of the filter is more complex than PDAF but the filter is 


again linear since INS error model is used. 


Batch size of the DTED area concemed can be changed. Both larger 
DTED areas for acguisition mode or smaller DTED areas for 


tracking modes can be selected using the same TSF structure. 


TSF gives solutions for various tracks selected. Actually, all tracks 
converge to the same index of the DTED grid (i.e. solution grid). 
However, for smooth terrains, there exist more than one position 
solution index and the tracks can be investigated separately in order 


to give more than one but finite number of navigation solutions. 


Results of the filter are good for both recursive and batch algorithms. 
The results compared with SITAN and TERCOM algorithms will be 


discussed in the simulations section. 


Simulation results of the implemented TSF will be presented at the 


“Simulations” section of the chapter. 


Simulations 


3.4.1. Simulation Model Development 


After the implementation of PDAF and TSF for TAN, the algorithms are 


tested using simple kinematic models. Kinematic models are prepared considering 
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the mid-course flight of a cruise missile with constant heading and velocity motion 


at constant altitude. 


Simulations for the TAN models are done using Simulink [58]. In order to 
perform simulations, first trajectory and INS kinematic models are formed. Then, 
DTED database model is prepared. Finally, SITAN, TERCOM, PDAF and TSF 
models are formed. Finally, the overall architecture is formed in order to perform 


simulations for position errors along east and north directions of the vehicle motion. 


In Figure 49, a general Simulink model for the studied TAN models is 
given. Loosely coupled integration structure is used for TAN models where INS is 
not updated at each TAN correction step but updated at a greater period. This is 
done in order not to influence INS results from possible fault corrected TAN 
solutions. Details of the simulation sub-models will be given in the following 


sections. 


INS MODEL 
TRAJECTORY 
MODEL 




















TERCOM ERROR 
STATES 
SITAN ERROR 
STATES 
PDAF ERROR 
STATES 
TSF ERROR 
STATES 
















DTED HEIGHT 
MODEL 


Measurements 
(Height Differences) 


o 
Ee 
= 
= 


INS ERROR = 
STATES 


Figure 49. General Simulink Model for TAN Models 





Navigation States 
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3.4.1.1. Trajectory and INS Model 


For the simulations, the motion of the vehicle is modeled considering the 
mid-course flight of a cruise missile with constant heading and velocity motion at 
constant altitude. Actually, this assumption is almost valid for mid-course flight of a 
generic cruise missile. For the INS model, white noise terms are added to velocity 
and position terms. Since the height terms will be taken from the DTED database 
according to vehicle’s latitude and longitude (i.e. horizontal positions), height is not 


considered in the vehicle’s state. 


Trajectory model of the vehicle considering continuous states can be 


modeled as follows: 
Krai (©) = F(X, (1) (3.43) 
where, 


Xi © =| Nya FE. 


traj? traj? 


T 
hy YN yaj3 VE yi VP ra | 


trap? VİN raja VE rrajo 
raj? Northward position of the vehicle 
rE, ,: Eastward position of the vehicle 
Altitude of the vehicle 


: Northward velocity of the vehicle 


vE: Eastward velocity of the vehicle 


traj 
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vD,,,: Down velocity of the vehicle 


traj * 


F(t)= 





= e O € o O 
= € O € € O 
= € O € c O 
= © O CO Sn 
= © O CO — O 
= © CO re O 





For the simulations, vehicle horizontal velocity and altitude are taken as 


constant values. These assumptions can be summarized as follows: 


VN „a; = Const. (3.44) 
VE ya = Const. (3.45) 
N; = Const. (3.46) 
VD ai =9 (3.47) 


In the same manner, for INS model white noise terms are added to the 


trajectory model: 


ae (1) = F(t) ys (t) + w(t) (3.48) 


where, 
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= T 
ys (0) = [FN ms PE ns Ams VN ms s VE mgs VP ms ] 


T 
WE) = wavs O; Wiep O; Wing (OX Wans (Os Warp (Oi KON 


w,(t)=N(0,07): Zero mean normal distribution with variance o; 


corresponding to related position and velocity 


For simulations with TAN algorithms, DTED are used in order to determine 
height differences at each time step as measurements which are given in eguation 
(3.11). DTED heights are given as a function of longitude and latitude. However, 
simulations are performed considering Cartesian coordinates in order to visualize 
navigation errors better in the simulations. Actually, latitudes and northward 
positions and longitudes and eastward positions are correlated. Correlations are 


defined as follows [73]: 








Sa ™ 3.49 
dt Ry +h 34) 
d vE 

H (3.50) 


dt" (R, +h)-cos4 

where, 

A: Latitude of the vehicle 
u: Longitude of the vehicle 


h: Height of the vehicle above ground 
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R,: Earth's polar radius 


R,: Earth’s equatorial radius 


The geometry of the earth is considered as an ellipsoid in WGS-84 


coordinate system. The earth’s polar radius R, and 


defined as [74]: 





ae a- (l-e) 
“ (-e -sin? 4)? 


2 + 2 1/2 
R, =a/(-e -sin” A) 
where, 
a: Semi-major axis 


e: First eccentricity squared 


WGS-84 values of these parameters are [74]: 


a =6378137+2m 


e” = 0.00669437999013 
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equatorial radius R, are 


(3.51) 


(3.52) 


(3.53) 


(3.54) 


R, and R, terms vary with changing latitude. However, in the simulations, 


changes in latitude will be less than one degree. Therefore, these terms are taken as 


constant values with their initial conditions taken in the simulations. 


Then, correlations for the white noise terms for longitudes and latitudes can 


be expressed as follows: 


ZE On 

A (Ry +1)? (3.55) 
2 

2 Or (3.56) 





Oo = 
“ [(R, +h)-cos Af 


In the simulations, discrete models are used considering the actual case. 
Height measurements will be taken at discrete intervals. Therefore, the models are 


discretized as follows: 
O(k) =1+ F(t)-T (3.57) 
where, 
D(k): State transition matrix 


T: Sample time (In the simulations, 7 =1 second taken.) 


Finally, trajectory model becomes: 
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Xing (k +1) = DA), (A) (3.58) 




















[7N,,(k+D] [1 0 0 T 0 0][7N,, (k)] 
rE, g(k+l)| |0 1 0 0 7 0||7E,,(k) 
h,y(k+D | (0 01 0 0 T|| hy 
| E . (3.59) 
VN, y(K+D)| |0 0 0 1 0 Of} | Wuk) 
vE„(k+1)| (0 0 0 0 1 0]]|vE, (k) 
| YD,y(k+D| [0 0 0 0 0 1]|vD,5(k)| 


Here, it should be noted that ®(k) term becomes constant considering 
constant velocity motion with the given sample time. Hence, linear discrete model 


is formed. In the same manner, for the INS model: 
Xis (k +1) = DA) Xs (1) + WA) (3.60) 
where, 


= T 
Xiys (kK) = [FN ins SPE ng Ayes VN wsi VE miss VP wns | 


WK) =| wy, E); We 


INS 


H); Wig ED Waring (E); Mot: D5 roy, OO] 


INS INS 


w(k)=N(0,07): Zero mean normal distribution with variance o? 


corresponding to related position and velocity 


White noise terms can be written as follows then: 
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2 
ON 


w, (k) = R Eye 


) 


2 
O 


m [(R, + iy: cos AP 





) 


Wy (k) = N(0, 07, ) 


W, (k) N0, 0) 


w, (k) = N(0,0;) 


Wy (k) = N(0, 055) 


Wye (k) = N(0,0;7) 


W,p(k) = N(0, 055) 


(3.61) 


(3.62) 


(3.63) 


(3.64) 


(3.65) 


(3.66) 


(3.67) 


(3.68) 


Here, it should be noted that, a, term is the variance of the barometric 


altimeter of the INS used. Radar altimeter errors will be added as measurement 


errors in the TAN algorithms. 


3.4.1.2. DTED Database Model 


As stated in the previous section, height of the vehicle is determined from 
the DTED maps according to the related latitude and longitude of the vehicle. For 
this purpose, Simulink “Lookup Table” blocks are used [58]. Then, measurement 
height differences are taken for the TERCOM, SITAN, PDAF and TSF models by 
adding INS white noises as system noise and radar white noises as measurement 


noise considering Figure 50. 





Figure 50. DTED Height Measurement Difference 


DTED height difference for the INS position was derived in eguation (3.11). 
Recall eguation (3.11) in order to investigate DTED model: 
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oh(k) = [pre (Lays (E), Ais O) + Wins | 


INS Height Model 


n [- raj (4), Araj (K))+ W adar (£) | 


Radar Height Measurement Model 








Here, it should be noted that first part of the eguation is the model of the INS 
height and second part is the radar height measurement model. Hence, the 
difference gives height difference term, dh(k). Using longitude and latitude values 
of the trajectory and INS models, height difference term can be obtained for 


simulations. 


Height difference, öh(k) parameter is used in SITAN model directly. 
However, batch height difference values are required for the applications of 
TERCOM, PDAF, TSF and “SITAN with Bank of Kalman Filters”. This is done by 
calculating all height difference terms in the 36 horizontal error bound of the INS 


and was derived in equation (3.12). Recall equation (3.12): 


oh,(k) = [pre 0,40) + Wys (k)| 


Height Model for All Positions in the DTED Batch 


[horeo lly Ary (BD) + War İİ 


Radar Height Measurement Model 








i=1,...,m(k) 


where; 


i: Index of the DTED grid node, 
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mk) 


Considering equation (3.12), formation of a batch of height differences for a 


: DTED grid size (selected as square of an odd number for INS 


position to be at the center of the DTED grid) 


5x5 DTED grid size is shown in Table 17. 


Table 17. Batch of Height Differences Formation for 5x5 DTED Grid 


A= +243" 
A= An +3" 
A = hs 
A= Ans 3" 
Ana, 23" 


Time:t=k M= His —2*3" M= ys —3" H iş H Hms t3" H= Mis + 2*3" 


Note: Index i=13 gives h(k) at position (Å ys » Hs ) Of INS, at time “t =1,”. 


Index i=1 gives dh(k) at position (4,,—2*3", 4s —2*3") of INS for DTED Level 


1, at time & . In the same manner, batch of öh(k) ’s are obtained. 


DTED database model formed in Simulink is shown in Figure 51. Here, 
latitudes and longitudes of the INS and trajectory model are determined first. Using 
*Lookup Table” blocks of Simulink, DTED heights for the INS and trajectory are 
found. Then height differences h(k) are determined in order to use in the TAN 
models. Moreover, for SITAN process, terrain slopes along longitude and latitude 


directions, (1,,h,) are also determined using the same procedure above. For 


“SITAN with Bank of Kalman Filters” model, batch of slopes are also formed 









































i=5 i=10 i=15 i= 20 i= 25 
i=4 i=9 i=14 i=19 i= 24 
i=3 i=8 i=18 i= 23 
i=2 i=7 i=12 i=17 i=22 
izl i=6 i=11 i=16 i=21 








considering batch formation explained in Table 17. 
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Figure 51. Simulink DTED Database Model 


For the simulations, DTED Level 1 data were reguired and they have been 


places of Turkey. 


obtained from HGK. The properties of DTED prepared for Turkey were given in 
Table 3 [16]. Horizontal accuracy of Level 1 DTED is defined as +130 m, and 
vertical accuracy as +30 m. On the other hand, horizontal accuracy of DTED Level 
2 data is +26 m. In the simulations, DTED Level 1 data are used considering the 
horizontal accuracies of DTED Level 2. Moreover, simulations with real DTED 


Level 2 data are also performed which were also obtained from HGK for a few 


In order to read DTED files for simulations, Matlab “Mapping Toolbox” is 


with the help of “Mapping Toolbox”. 


used [75]. Hence, binary DTED files are directly used as text files in the simulations 
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3.4.1.3. TERCOM Model 


TERCOM was studied in detail in Chapter 2 with the simulations performed 
and the results were discussed. For TAN algorithm simulations, the models for 
TERCOM are formed considering the well-known MAD and MSD processes 
derived for TERCOM. 


Recall eguations derived for MAD and MSD processes given from eguation 
(2.9) to (2.11): 


MSD, =(1/ NY3.(S, Sy) 
izl 


N 
MAD, = VM), 





Sy- Sa 


PIC <C il, where a minimum of C, is sought, 
Sa = P[C, >C], where a maximum of C „ is sought. 


Examination of the expressions for the MAD and MSD processors indicates 
that both of these correlators can be viewed as distance measures, where the 
dimensions of the space for which these distances are defined correspond to the 


number of elements in the profiles. 


In order to form TERCOM model in Simulink, “Signal Processing Blockset” 
is used [76]. In this blockset, RMS, mean and minimum selection operations 


defined for MAD and MSD processes can be done in real-time. 
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TERCOM model formed in Simulink is shown in Figure 52. Here, at the 
upper part MAD process, at the lower part MSD process occurs. Using batch of 
height differences “d6h(k) “obtained from DTED database model, minimum of 
MAD and MSD functions are determined. It should be noted that, TERCOM 
Simulink model works in real-time by calculating MAD and MSD functions at each 
time step. In actual applications, these functions are calculated once after some time 
of the operation begins. However, considering calculated MAD and MSD functions 


at some definite times of operation, TERCOM processes can be understood. 


On the other hand, position corrections are done considering calculated 
position indices from TERCOM. As a result of this, position accuracies for 
TERCOM are within the limit of the DTED grid size. In other words, TERCOM 


horizontal position accuracies can not be better than the DTED grid accuracy used. 











MATLAB 
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del DTED 
MATLAB Fen 


Batch 
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Mn. MAD (Yellow) Correction 
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Running 
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MSD 
Correction 






Mnmum1 


MAD lat nd.(yeliow) 
MAD jon nd. (magenta) 


Figure 52. Simulink TERCOM Model 
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3.4.1.4. SITAN Model 


SITAN was studied in detail in Chapter 2 with the simulations performed 
and the results were discussed. For TAN algorithm simulations, the models formed 


for SITAN are directly used. Two models are used for SITAN: 


1. Standard EKF for SITAN for tracking mode (Errors less than 100 
meters for DTED Level 1 maps) 


2. Bank of EKF’s for acquisition mode (Errors grater than 100 meters 
for DTFD Level 1 maps) 


For tracking mode, eguations derived for SITAN given in Table 10 are used. 
EKF eguations are written in Simulink using S-functions [58]. An S-function is a 
computer language description of a Simulink block. S-functions use a special 
calling syntax that enables the user to interact with Simulink eguation solvers. This 
interaction is very similar to the interaction that takes place between the solvers and 
built-in Simulink blocks. The form of an S-function is very general and can 


accommodate continuous, discrete, and hybrid systems [14]. 


Simulink model with standard EKF for SITAN for tracking mode is shown 
in Figure 53. Here, EKF uses the height difference between the INS and the 
trajectory and the slopes at the INS position. 


For acguisition mode, parallel EKF structure with 3x3 and 5x5 grid size is 
used. In order to select position fix, eguation (2.40) is used. Recall eguation (2.40) 


which gives AWRS value of the selected filter: 
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INS lation 
Error 


INS Error 
States. 


Param's 
SITAN Error 
States 


SITAN EKF 


Figure 53. Simulink Single EKF SITAN Model 





1/4 A, 
AWRS ju fiter = N 27 PH +R, 


jth filter 


where; 


AWRS inttr: Average Weighted Residual Squared of the j’th filter, 


This AWRS value is the average weighted residual squared between the 
predicted ground clearance for each filter and the ground clearance measured by the 
radar altimeter for each time ¢,. By examining the minimum AWRS values for each 
filter after a sufficiently large number of measurements have been processed, the 


correct filter and its associated state error estimates are chosen. 


Formation of the SITAN model for acquisition mode is quite complex. In 


Figure 54 and Figure 55, SITAN models with bank of EKF’s are shown. 
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KFSS 
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Bank of SITAN INS States 
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KF Selection 





Bank of KFs 


KF Index 


Min. AWRS 


Figure 54. Simulink Bank of EKF SITAN Model for 5x5 Grids 


In these models, parallel Kalman filters run with different initial conditions 
and slopes. Single EKF model is used for tracking mode simulations, where bank of 


EKF’s models are used in the acquisition mode simulations. Simulation results of 


SITAN filters for tracking and acquisition modes were also given in Chapter 2. 


176 





























































































































































































































































































































INS3İ a KF States INS32 z 
7 F 

EKFbank iel EKFbank (ef: 

: KF32 7 

































































Val I 
—— 
Min. AMRS lax 


Minmum 


KF Index 


Figure 55. Simulink Bank of EKF SITAN Model for 3x3 Grids 








3.4.1.5. PDAF Model 


PDAF model is formed in Simulink considering the equations in section 
3.2.2.2 which were summarized in Table 15. Here, it should be noted that Simulink 
architecture is very simple. Only a single PDAF is used in all of the simulations. 
Batch size of the filter (i.e. considered grid size) can be changed independent of the 
PDA filter. PDAF model can be used for both acquisition and tracking modes of the 
TAN solution. 


In the PDAF model, following processes are performed: 
1. Height differences are taken from the DTED database model. 


2. Gating process is done in order to extract impossible position 


solutions. 
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3.4.1.6. 


3. Height differences are averaged in order to form past measurement 


information. 


4. Averaged height differences are used in the PDAF in order to 


determine PDAF error states. 


Simulink PDAF model is shown in Figure 56. 


MATLAB Delay I i 
Line a jeu <2 
del DTEDs PDAF 


Measurement PDAF ie 
Gatng Past Measurement Mean Unbuffer 
Information 


Figure 56. Simulink PDAF Model 


TSF Model 


TSF model is formed in Simulink considering the eguations in section 3.3.3 


which were summarized in Table 16. Here, it should be noted that again Simulink 


architecture is very simple like PDAF model. However, more complex operations 


are done in TSF model, since bank of Kalman filter operations are performed. Bank 


of Kalman filter operations are performed using an S-function. Batch size of the 


filter (i.e. considered grid size) can be changed independent of the filter. TSF model 


can be used for both acguisition and tracking modes of the TAN solution as 


explained in the related sections before. 
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In the TSF model, following processes are performed: 


Height differences are taken from the DTED database model. 


Gating process is done in order to extract the impossible position 


solutions. 


Height differences are used in the TSF in order to determine TSF 


error states. 


Simulink TSF model is shown in Figure 57. 





del DTEDs TSF Sistes 


Messurement TS Filter 
Gsting 


Figure 57. Simulink TSF Model 


3.4.2. Case Studies 


Case studies with simple kinematic models are performed for three different 


cases: 


1. Simulations with DTED Level 1 


2. Simulations with DTED Level 2 
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3. Simulations with Various DTED Grid Sizes 


For various DTED level simulations, different terrain types are selected. 
Then, simulations for tracking and acguisition modes are performed simultaneously 


using the implemented TAN algorithms compared with the well-known algorithms. 


3.4.2.1. Simulations With DTED Level 1 


3.4.2.1.1. Terrain Selection 


Simulations are performed for three different types of terrains: 


1. Rough terrain 


2. Smooth terrain 


3. Mountainous terrain 


In order to determine the reguired terrains, Microdem/ TerraBase II 
Software [77] is used. Terrain contours of 50 meters and the trajectory paths 
obtained from the software is shown in Figure 58. Then, terrain heights versus time 
plots for the selected terrains are given in Figure 59, Figure 60 and Figure 61. 
Finally, terrain parameters are calculated for the selected terrains and summarized 


in Table 18. 
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Figure 58. Terrain Contours for TAN Simulations 
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Figure 59. Terrain Height vs. Time for Rough Terrain 
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Figure 60. Terrain Height vs. Time for Smooth Terrain 
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Figure 61. Terrain Height vs. Time for Mountainous Terrain 


Table 18. Terrain Parameters for TAN Simulations 




















Terrain Type Rough Smooth Mountainous 
Mean height of the terrain profile 1093 m 1104 m 1177 m 
Sigma-T 71.9 m 34.1 m 212.9 m 
Sigma-Z 162 m 3.7 m 23.1 m 

X, 674.6 1309 m 1302 m 
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3.4.2.1.2. INS Model Verification 


In order to perform simulations, verification of the INS model used is 
reguired first. For cruise missiles, generally 1.0 nm/hr INS guality is reguired. In the 
simulations, simple INS model is formed considering eguation (3.60). In this 
eguation, white noise terms are added to positions and velocities considering 


constant velocity flight. 


In Table 19, simulation parameters used for the INS model are given. Using 
the values provided, INS model is tested for horizontal position errors. Horizontal 


position and velocity errors are given in Figure 62 and Figure 63. 


Table 19. INS Model Parameters for 1.0 nm/hr Ouality 

















Initial vehicle velocity 240 m/s 
Initial INS east velocity bias 0.5 m/s 
Initial INS north velocity bias 0.5 m/s 
INS horizontal position standard deviation (0,,, 0.) 9m 
INS velocity standard deviation (oy, 0,;) 0.05 m/s 























As it can be seen from the horizontal position errors, INS quality is about 3.0 
nm/hr. Actually, selecting worse INS guality than the real system used is preferred 
in order to test the performance of the TAN algorithms used. In the simulations, 
small operation times (like 100 seconds) are used for TAN applications where INS 


is not updated during the simulation period. 
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Horizontal INS Position Errors vs. Time 
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Figure 62. 


Horizontal INS Velocity Errors vs. Time 
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Figure 63. Horizontal Velocity Errors of the INS Model Used 
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3.4.2.1.3. Simulation Results 


Simulations are performed for two modes of operation of TAN algorithms: 


1. Tracking mode, where SITAN single filter, PDAF and TSF are 
compared with Monte Carlo simulations along the trajectory 


(Recursive Solution); 


2. Acguisition mode, where PDAF, TSF and TERCOM are compared 
with Monte Carlo simulations for the position update at a defined 


time (Batch Solution). 


3.4.2.1.3.1. Simulations for Tracking Mode 


First, simulations for tracking mode are done. Parameters used in the 
simulations are given in Table 20. Monte Carlo simulations of 100 runs are 


performed and the following plots are obtained for SITAN, PDAF and TSF. 


1. Northward and eastward position errors; 


2. RMS values of north and east positions. 


For tracking mode simulations, TERCOM is not used. In order to apply 
TERCOM algorithm, larger DTED grid size and large initial position errors are 
reguired. When TERCOM algorithm is applied for small grid size, false position 
fixes occur with high percentages since INS error model is not used in TERCOM 


algorithm. 


185 


Table 20. Simulation Parameters for Tracking Mode 






































Initial INS position deviation (one axis) 60 m 
Initial vehicle velocity 240 m/s 
Initial INS east velocity bias 0.5 m/s 
Initial INS north velocity bias 0.5 m/s 
INS horizontal position standard deviation 9m 
INS altitude position standard deviation 3m 
Radar altimeter standard deviation 3m 
INS velocity standard deviation 0.05 m/s 
DTED Grid Size (for PDAF and TSF) 3x3 











Terrain Type 1 (Rough Terrain): 


Delta rN vs. Time 














Figure 64. Northward Position Error vs. Time 
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Figure 65. Eastward Position Error vs. Time 
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Figure 66. Northward Position RMS Error vs. Time 


187 


RMS Delta rE vs. Time 
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Figure 67. Eastward Position RMS Error vs. Time 
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Figure 68. Total Position RMS Error vs. Time 
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Figure 69. Northward Position Error vs. Time 
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Figure 70. Eastward Position Error vs. Time 
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Figure 71. Total Position RMS Error vs. Time 
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Figure 72. Northward Position Error vs. Time 
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Figure 73. Eastward Position Error vs. Time 
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Figure 74. Northward Position RMS Error vs. Time 
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RMS Delta rE vs. Time 
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Figure 75. Eastward Position RMS Error vs. Time 
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Figure 76. Total Position RMS Error vs. Time 
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Simulations for tracking mode are done for a small period of time (i.e. 100 
seconds for tracking mode) in order to visualize the performance of the 
implemented TAN algorithms. In the actual navigation system of a cruise missile, 
INS will be updated at discrete time intervals according to the TAN algorithm used. 
In the simulations performed above, after 10 to 40 seconds of operation, since INS 
is not updated errors grow and TAN algorithm can not be used since small DTED 
grid size is selected. TAN performance depends on the terrain type. For 
mountainous terrains, implemented TAN algorithms find position fixes faster than 
rough terrains. Hence, INS update time can be determined from the selected terrain 


properties. 


From the simulations, it is seen that better results than SITAN are obtained 
for rough and mountainous terrain types. For the smooth terrain, SITAN seems to 
show better results. However, response of the SITAN filter is also not stable and 
navigation solution cannot be obtained for smooth terrain. TSF and PDA filter 
results are considerably good, since if navigation solution does not exist, the filters 
follow INS error model which is actually a desired feature. From the Monte Carlo 
simulations, position RMS errors of the TSF and PDAF algorithms become less 
than 50 meters for mountainous terrains; in other words, a decreased navigation 
error is obtained. As it can be seen from the simulation results TSF behaves as a 


correction shift along the INS error model. 


3.4.2.1.3.2. Simulations for Acguisition Mode 


Parameters used in the simulations for acguisition mode are given in Table 
21. Monte Carlo simulations of 100 runs are performed and northward and eastward 
position errors at the update time are obtained from the plots for TERCOM, PDAF 
and TSF. 
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Table 21. Simulation Parameters for Acguisition Mode 






































Initial INS position deviation (one axis) 400 m 
Initial vehicle velocity 240 m/s 
Initial INS east velocity bias 0.5 m/s 
Initial INS north velocity bias 0.5 m/s 
INS horizontal position standard deviation 9m 
INS altitude position standard deviation 3m 
Radar altimeter standard deviation 3m 
INS velocity standard deviation 0.05 m/s 
DTED Grid Size (for PDAF and TSF) 11x11 











Terrain Type 1 (Rough Terrain): 
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Figure 77. Northward Position Error vs. Time 
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Figure 78. Eastward Position Error vs. Time 
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Figure 79. TSF Indices vs. Time 
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Terrain Type 2 (Smooth Terrain): 
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Figure 80. Northward Position Error vs. Time 
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Figure 81. Eastward Position Error vs. Time 
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Terrain Type 3 (Mountainous Terrain): 
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Figure 82. Northward Position Error vs. Time 


Delta rE vs. Time 


> J 





Figure 83. Eastward Position Error vs. Time 
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TSF Indeces vs. Time 














Figure 84. TSF Indices vs. Time 


From the simulations, it is seen that similar results with TERCOM are 
obtained for rough terrain type. For the smooth terrain, both algorithms do not have 
navigation solution. As it can be seen from the results, TSF and PDA filter can be 
used also for acguisition mode of TAN solution using the considered DTED size. 
SITAN bank of Kalman filters is not considered for acguisition mode here; since, 
for the initial error given 121 Kalman filters should be run for the simulation. On 


the other hand, using TSF and PDAF, same solutions with TERCOM are obtained. 


In the simulations, TSF indices for navigation solutions are also given. Here, 
best definite number of navigation solutions converges to the same position index in 
a few seconds especially for mountainous terrain. However, for smooth terrain, 
indices changes in time unlike other terrain types. On the other hand, pruning 
method is selected such that best definite number of tracks with minimum 
likelihood functions selected. Working on other pruning methods may improve 


navigation solutions for smooth terrains. 
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Another critical point in the simulations is the percentage of false fix for 
acguisition mode. In TERCOM, since INS error model is not considered for 
correlation process, there is always a probability of false fix in the position 
solutions. In the Monte Carlo simulations performed, a few false position fixes 
occurred for TERCOM for rough terrain type. On the other hand, with PDAF and 


TSF no false position fixes occurred. 


3.4.2.2. Simulations With DTED Level 2 


3.4.2.2.1. Terrain Properties 


Simulations with DTED Level 2 are performed for a single terrain type 
shown in Figure 85. Parameters of the selected terrain showed the area to be a rough 
terrain. Then, terrain heights versus time plot for the selected terrain is given in 


Figure 86. 





Figure 85. DTED Level 2 Terrain for TAN Simulations 
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Figure 86. Terrain Height vs. Time for DTED Level 2 Terrain 


3.4.2.2.2. Simulation Results 


Simulations are performed for both tracking and acquisition modes of 


operation as in DTED Level 1 simulations. 


3.4.2.2.2.1. Simulations for Tracking Mode 


Monte Carlo simulations of 100 runs are performed and position errors are 
obtained for SITAN, PDAF and TSF. Parameters used in the DTED Level 2 
simulations are given in Table 22. It should be noted that since DTED Level 2 is 


used, initial INS error is selected smaller considering grid size of the DTED used. 


Table 22. Simulation Parameters for DTED Level 2 Tracking Mode 








Initial INS position deviation (one axis) 25m 





Other parameters Same as in Table 20. 
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Delta rN vs. Time 
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Figure 87. Northward Position Error vs. Time 
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Figure 88. Eastward Position Error vs. Time 
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RMS Delta rN vs. Time 


i 
i 
i 
i 
i 
D 
D 
e 


50 þ----------- 


- — : SITAN 





Figure 89. Northward Position RMS Error vs. Time 
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Figure 90. Eastward Position RMS Error vs. Time 
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RMS Total Position Error vs. Time 
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Figure 91. Total Position RMS Error vs. Time 


From the simulations, it is seen that better results than SITAN are obtained 
for the selected terrain type. From the Monte Carlo simulations, position RMS 
errors of the TSF and PDAF algorithms become less than 25 meters for the selected 
terrain. It should be noted that DTED Level 2 grid size accuracy is about 30 meters. 


Hence, real-time accuracy is increased. 


3.4.2.2.2.2. Simulations for Acguisition Mode 


Parameters used in the simulations for acguisition mode are given in Table 
23. Monte Carlo simulations of 100 runs are performed and northward and eastward 
position errors at the update time are obtained from the plots for TERCOM, PDAF 
and TSF. 
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Table 23. Simulation Parameters for DTED Level 2 Acguisition Mode 








Initial INS position deviation (one axis) 140 m 











Other parameters Same as in Table 21. 
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Figure 92. Northward Position Error vs. Time 
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Figure 93. Eastward Position Error vs. Time 
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TSF Indeces vs. Time 











Figure 94. TSF Indices vs. Time 


Simulation results show that TAN algorithms also work with DTED Level 2 
for acguisition mode. However, DTED Level 1 results seem to be better than DTED 
Level 2 results. Unfortunately, there were not sufficient DTED Level 2 maps for 
simulations in order to compare simulation results in detail. Actually, vehicle 
velocity directly influences TAN performance. For cruise missiles, DTED Level 1 
maps are sufficient for mid-course flight navigation solution where INS position 
fixes less than 50 meters can be obtained. For faster vehicles like cruise missiles, 
rapid changes in the terrain profile as in DTED Level 2 decreases TAN 
performance. As a result of this, use of DTED Level 1 maps for TAN acguisition 


mode seems to perform better solutions. 
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3.4.2.3. Simulations with Various DTED Grid Sizes 


Final part of the case studies is done with various DTED grid sizes for 
PDAF and TSF. Here same initial position errors are taken for simulations in 
tracking mode along rough terrain. Parameters used in the simulations for 
acquisition mode are given in Table 24. Monte Carlo simulations of 100 runs are 


performed and northward and eastward position errors are obtained for PDAF and 
TSF. 


Table 24. Simulation Parameters for Various DTED Grid Sizes 








Initial INS position deviation (one axis) 80 m 
DTED Grid Size (for PDAF and TSF) 3x3 


5x5 
7x7 
9x9 
11x11 
Other parameters Same as in Table 21. 





























PDAF: Delta rN vs. Time 
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Grid Sizez 7 
Grid Size= 9 
Grid Size= 11 
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Figure 95. PDAF Northward Position Error vs. Time 
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It should be noted that initial position errors are taken small in order to have 
solutions with small DTED grid sizes. Hence, effects of selecting larger DTED grid 


sizes are examined in this section. 


PDAF: Delta rE vs. Time 
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Figure 96. PDAF Fastward Position Error vs. Time 


PDAF: RMS Total Position Error vs. Time 
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Figure 97. PDAF Total Position RMS Error vs. Time 
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TSF: Delta rN vs. Time 
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Figure 100. TSF Total Position RMS Error vs. Time 
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It is seen that changing grid size for the same initial position errors for TSF 
slightly changes simulation results. Again errors are bounded and limited with the 
related grid solution. However, selecting larger grid sizes for PDAF solutions 
generally increase position errors. This is due to PDA procedure where weighted 
averages of the all grid points are taken into account for navigation solution. 
Therefore, it can be concluded that PDAF DTED grid size should be selected in 


accordance with the position errors. 


3.4.2.4. Discussion 


From the simulations performed, several conclusions are achieved about the 
implemented TAN algorithms. The advantages of the new algorithms proposed can 


be summarized as follows: 


1. Real-time TAN solution can be obtained with a single PDA filter. 
Since past measurements are taken into account, by changing the 
buffer size of the measurements the filter, measurements are 


smoothed. 


2. Real-time TAN solution can be obtained with a single TSF structure. 
However, TSF operations are more complex than SITAN. On the 
other hand, in TSF, more than one track is selected in order to 
determine navigation solution. Hence, probability of false fix 


decreases unlike TERCOM. 


3. Real-time TAN solution is obtained by considering horizontal 
position errors of DTED used in real-time PDA filter and TSF. 
Hence, horizontal position states are added to the Kalman filters used 


in PDAF and TSF. 
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. Application of the filters is simple and the filters are linear, since 


INS error model is used. 


Batch size of the DTED area concerned can be changed independent 
of the model used. Both larger DTED areas for acguisition mode or 
smaller DTED areas for tracking modes can be selected using the 


same filters. 


Results of the filters are good for both recursive and batch 
algorithms. For tracking mode, position RMS error is less than 50 
meters. Moreover, PDAF shows stable response. For smooth terrains 
where no navigation solution exists, PDAF follows the INS error 


model which is actually a desired feature. 


TSF can be considered as a real-time TERCOM process for large 
position errors, i.e. large DTED batch size. Possibility of false 
position fixes decrease with TSF when compared with TERCOM. 
On the other hand, for small position errors, decreasing the 
weighting factor of the past measurements for TSF, better real-time 


solutions can be obtained. 
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CHAPTER 4 


CASE STUDY 


In this chapter, simulation results of the implemented TAN algorithms are 
presented for a cruise missile model. First, a 6 DOF simulation tool is developed in 
order to model cruise missile mid-course flight. Then, related sub-systems are 
modeled in order to reflect mid-course flight controls and cruise missile navigation 
system error models. Then, simulations are performed with PDAF and TSF TAN 
models with actual flight conditions. Finally, simulation results are compared with 
major TAN algorithms considering other flight parameters of the cruise missile 


model. 


4.1. Simulation Tool Development 


The simulation tool developed for the cruise missile is capable of 
performing full mid-course flight simulation of the cruise missile modeled. 
Actually, a generic simulation tool applicable to all air vehicles is considered except 


for guidance methods applied. 


In order to investigate the performance of the TAN algorithms improved, a 


realistic 6 DOF simulation tool is reguired. 6 DOF simulation model is built in 
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Simulink [58] and mid-course flight of the cruise missile is simulated with the 
model developed. General 6 DOF simulation model architecture is shown in Figure 
101 where TAN algorithms are used with the loosely coupled architecture for aiding 
INS. 





AIRFRAME 






AUTOPILOTS ACTUATORS 
INS Accelarations and 


Angular Velocities MISSILE 
SENSORS 










TAN MODELS 
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Figure 101. General 6 DOF Simulink Model with Implemented TAN Models 


6 DOF simulation model is formed from the following main sub-systems: 
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1. Airframe: 


2. Autopilots: 


3. Actuators: 


4. Sensors: 


5. Guidance: 


Airframe is composed of various models in order to 
simulate the dynamic behavior of the cruise missile. 
Dynamics, kinematics, aerodynamics, propulsion, force 
and environmental modules constitute airframe sub- 


system. 


A variety of controllers are implemented for mid-course 
flight of the cruise missile: roll control, pitch 
acceleration controller for altitude hold, yaw stability 
augmentation, and BTT (bank-to-turn) heading angle 
tracker autopilots. All autopilots are derived by the pole 


placement technigues. 


Second-order actuators with rate and position limiters 
are used in the 6 DOF model in order to control the 


elevator, rudder and aileron. 


Strapdown INS sensor errors are modeled in order to 
reflect bias, drift, scale factor and misalignment errors 
for accelerometer and gyro outputs. Moreover, 
barometric and radar altimeter outputs are also 


modeled. 


For the selected waypoints along the missile path, 
heading correction is applied using the heading angle 


tracker. 


6. INS Error Model: o Strapdown INS error model given in equation 


(2.6) is used in the simulations in order to 
reflect INS velocity and position errors 


considering improved bias and drift models of 
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the INS sensors. However, altitude channel of 
the INS model is not used considering 


barometric altimeter measurements. 


7. TAN Models: Derived TAN models in the previous chapters for 
TERCOM, SITAN, PDAF and TSF are directly used 
in the 6 DOF simulation model by tuning parameters 


of the derived filters. 


6 DOF simulation model sub-systems will be discussed in the following 
sections except for INS error and TAN models which were investigated in detail in 


the previous chapters. 


4.1.1. Airframe 


Airframe sub-system models dynamic behavior of the cruise missile model. 
6 DOF eguations of motion are derived and used in this section with the related sub- 


systems described in the previous section. 


Newton’s law is applied for translational motion and Euler’s law is applied 
for rotational motion in order to model cruise missile dynamics over an elliptical 
earth considering earth's rotation. First, reference frames are defined. Then, 
Newton’s and Euler’s equations are derived considering forces and moments on the 


system. Finally, kinematic eguations are derived. 
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4.1.1.1. Reference Frames 


Coordinate frames are reguired in order to define the motion of the vehicle 
which is considered. Moreover, for kinematic eguations, they have to be defined 
carefully. When the rotation of the earth and earth's geometry is considered, various 
coordinate frames have to be defined. Coordinate frames used for the simulations 


can be classified as follows: 


1. Inertial Frame (Geocentric Inertial-J2000 Frame), 3, 


2. Earth Centered Earth Fixed (ECEF) Frame, 3, 


3. Geographic (North, East, Down — NED) Frame, 3, 


4. Body Frame, 3, 


5. Wind Coordinate Frame, 3, 


The position and attitude of the missile with respect to inertial frame is 
found using kinematic eguations. Since time of flight for the cruise missiles are 
guite long, the effects of carth's curvature and earth's angular velocity should be 


added in the transformations. 


In Figure 102 and Figure 103, these frames were presented. Transformation 
matrices are derived between these frames and using translational and rotational 
transformation eguations, positions and attitudes are defined in the proper reference 


frames as follows considering tensor algebra (78): 
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Figure 102. Inertial and ECEF Reference Frames 
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Figure 103. Geographic and Body Reference Frames 
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pla) — Êb) 70) (4.1) 


Gan = GEN GO 42) 


where, 


r: Position vector 
C: Transformation matrix (Direction Cosine Matrix, DCM) 
O: Skew symmetric matrix form of angular velocity vector, © 


ON: Angular velocity vector of the body frame with respect to inertial 


frame in body frame 


a,b,i: Indices of the related reference frames 


Then, transformation procedures for the defined coordinate frames can be 


applied considering “Rotated Frame Based” (RFB) seguences [78]: 


a. Transformation from Body to Geographic Frame 


= (2) = (m) z (n) 
ui ü ui 
S. 3 > $ 2 > 5, l > 3, 
v 0 b 


4.3 
cy:cO cy:-sO-sP-sy-cp cv:-sO-cp+swv:-sp ee) 


CE» =lsy-cO sy sO-sp+cy-cp sy-s0-ch—cy-sh 
— s0 c0- so c0- co 
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(v: Yaw angle, 0: Pitch angle, g: Roll angle, c: Cos, s: Sin) 


b. Transformation from Geographic to ECEF Frame 








o. >, > 3, 
KOS arty 
2 
N —sh-cd —su —cA-cu (4.4) 
CE9 -|-s4.su cu —cA-su 
cA 0 -sÀ 


( u: Longitude, 2: Geocentric Latitude) 


c. Transformation from ECEF to Inertial Frame 





cO-t -Qt 0 
CEP =| sQ-t Qt 0 
0 0 1 


(4.5) 


(© : Earth’s angular velocity) 


d. Transformation from Body to Wind Frame 
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ca:cB sp sa-cp 
CD =|-ca-sB cB -sa-sp 
-sa 0 ca (4.6) 


(a: Angle of attack, 4: Side slip angle) 


Finally, between different coordinate frames successive transformations can 


be done as follows [78]: 
CD = EH) Alen), ED (4.7) 


Ce) — Ölen Gem) — Ales), G0.07 (4.8) 


4.1.1.2. 6DOF Eguations of Motion 


For translational motion, Newton's law is applied with respect to inertial 


frame considering the conventions in Figure 103 in vector notation. 


m-D5j0=F,,+m-8 (4.9) 
where, 

m: Mass ofthe vehicle 

D,: Differential operator in inertial frame 

Viio: Velocity vector of the center of mass with respect to inertial frame 
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E, p: Aerodynamic and propulsion forces acting on the vehicle 


g: Gravitational acceleration 


Then, the equations are derived in geographic frame from equation (4.9) 


considering kinematic equations as follows: 


D 110 = DV 110 + Oai *V 110 (4.10) 


D, (@,,; x Fio) = DO, Xb 49 + Ooi X D Fio 
N 
ö,; ;=Const. 
O, X DF 110 = Oni *| Difyr t Diyo |= Oni *Vare 


Vaio = Vare Y Oeri *V410 


DV 110 = DW gig + 2+ Oai X Vage +O Daş XF 110 (4.11) 
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After replacing eguation (4.11) in (4.10), translational eguations of motion 


can be written in geographic frame. 


pt 82: Oi, XV aj Oui X Oe), XF110 (4.12) 


A 


Equation (4.12) is written in tensor form in geographic frame because of its 


use in scalar manipulations. 


1 
52) (gb) FO 4 g2) (2) tg) — E VE 
Viz = ire Ce i E +g“ = 2: Oyj, X Oli; * Oj, T, i (4.13) 


where, 


Ce) = Ce) COT: DCM from body to geographic frame 


FA : Aerodynamic and propulsion forces acting on the vehicle 
g = = [0 0 el: Gravitational acceleration vector 
m vE [Q-cosA 0 -Q-sin A]: Earth’s angular velocity in 3 


ve [vN vE vD] : Velocity of the vehicle in geographic frame 


Next, Euler’s law is applied for rotational motion about the center of mass of 


the vehicle. 
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DH ,=D,H,+@,,xH,=M, (4.14) 


where, 


H ,: Angular momentum vector about center of mass (H, =/,-@,,,) 


M ,: Moment vector about center of mass 


Again equation (4.14) is written in tensor form in body frame in order to 


obtain angular velocities of the vehicle. 


Tb) 0) (b) F) =) ayı) 
I, Oy, tOn de O =M, 


A =I (ah İP a) sa) (4.15) 
where, 

1 0 0 
I PSO Z 0 Moment of inertia matrix of the vehicle 

0 0 7 


MO): Aerodynamic moments acting on the vehicle in body frame 


o) = [ p g aes Roll, pitch and yaw rate of the vehicle 
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4.1.1.3. Kinematic Equations 


Translational and angular velocities of the vehicle are obtained from 6 DOF 
equations of motion. Then, positions and attitudes are derived using kinematic 
equations. 

Tae Diye = Vale (4.16) 


where, 


Fg >İPN rE rDJ: Northward, eastward and downward positions 


A (ib) Alib) ~(b) 
CC, 


ÖL) = G6DT Öle) (4.17) 


Using the elements of C8) matrix given in equation (4.3) roll, pitch and 


yaw angles (i.e. J, O and y ) are obtained. 


4.1.1.4. Aerodynamics and Propulsion 


In the previous sections, dynamic and kinematic eguations of motion are 
derived. In order to solve dynamic eguations of motion, forces and moments acting 


on the body should be known. External forces and moments acting on the body are 
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due to aerodynamics and propulsion. Forces and moments acting on the body are 


summarized as follows [79]: 


Cy Ja 
FO = FO +F” =7-S. C, |+| 0 (4.18) 
Cy 0 
= € p, 
FO = COM. FO + FM = COM) Z.S C, ||+] 0 (4.19) 
=C; 0 
C,:b 
M” =g-S-|C, e (4.20) 
CD 
where, 
FP: Aerodynamic forces acting on the vehicle in body frame 
FW: Aerodynamic forces acting on the vehicle in wind (stability) frame 
F = : — Propulsion force acting on the vehicle, (F, : Scalar value) 
; 1 
g: Dynamic pressure (7 = ee y?) 
p: Density of the ambient atmosphere 
V: Velocity of the vehicle 
S: Reference area (Theoretical wing area for aircraft type vehicle) 
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b: Span 


C: Chord 


C,,C,,C,: Aerodynamic force coefficients (Axial, side and normal force 


coefficients) 


C,,C,,C,: Aerodynamic force coefficients (Drag, side and lift force 
coefficients) 
C,C,,,C,: Aerodynamic moment coefficients (Rolling, pitching and 


yawing moment coefficients) 


Aerodynamic force and moment coefficients are modeled generally by 
simple Taylor series expansion, including only the linear terms and making all 
derivatives tabular functions of Mach number and angle of attack considering small 


side slip angles as follows [79]: 


C, =C, (M,a)+C, (M,a)-a+C, (M,a)-öe+ 

c (4.21) 
C M,a ge —— 
zi )-q Ww 


Cp =Cp, (M,4)+C, (M,a): a (4.22) 


C, =C;,(M,a)+C, (M,a): B+C, (M,a)-da+ 


rng ee Grip a, 2 
Ys, ? e , p 2y Y, , 57 
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C,=C,(M,a)+C,,(M, a): B+C, (M.a)-Sa+ 


C, (M,a)-6r+C (Ma): a G (M yi oon 
a A aa 2V 
Cn = Cy (M,A)+C,, (M, A): +C, (M,a)-de+ 
c (4.25) 
C M,a -g:-— 
m (M,&)-g ay 
C, =C, (M,a)+C, (M,a)-B+C, (M,a)-ör+ 
b (4.26) 


b 
C (M,a)-64+C (M,a): p:-— +C (M,a)-r-— 
„M,a):5a +C, (M,a)-p-5—+C, (M,a): r> 


The equations above represent a linear model for constant Mach numbers 
and angle of attack values. In force coefficients, the effect of body rates on the lift 


and side forces, i.e., the C, , Cy , and C, derivatives can often be neglected [79]. 


For moment coefficients, rolling and yawing moment coefficients have negligible 


trim coefficients C, and C, [79]. 


For 6 DOF simulations, aerodynamic coefficients are obtained from USAF 
Digital DATCOM software [80]. First, a solid model of a generic cruise missile is 
formed using a CAD software. Next, mass and inertia properties are obtained from 
the solid model. Finally, USAF Digital DATCOM model of the cruise missile is 
prepared in order to obtain aerodynamic coefficients and coefficient derivatives. In 
Figure 104, a view of the cruise missile solid model is shown. In Table 25, some 


properties of the cruise missile used in the simulations are given. 


226 





Figure 104. Cruise Missile Solid Model 


Table 25. Cruise Missile Model Specifications 












































Length: 6.25 meters 

Initial Missile Weight: 1200 kg 

Fuel Weight: 450 kg 

Diameter: 53.34 cm 

Wing Span: 4.64 meters 

Wing Chord: 0.689 meters 

Speed: 0.4 — 0.8 Mach 

Cruise Speed: 0.7 Mach 

Engine Type: 600 Ib-f (~2500 N) Type 








Some of the coefficients reguired in eguations (4.21) and (4.26) are not 
directly obtained from Digital DATCOM. Therefore, using the coefficients obtained 
from the Digital DATCOM manual [80], reguired coefficients are derived. 
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Propulsion force acting on the missile is considered from turbojet and acting 
along axial axis only. Turbojet force is obtained as tabulated results with respect to 
Mach number and altitude above sea level. The specific fuel consumption (SFC) b, 
is an important indicator for the efficiency of the turbojet. It is defined by the ratio 


of fuel flow to thrust as [79]: 


b, = = (4.27) 
p 

where, 

Mp: Fuel flow rate 


Required thrust F, is obtained from Mach control loop. Using look up 


tables, SFC and fuel flow rates are obtained. Hence, mass of the missile is modeled 


for simulations. 


4.1.1.5. Environmental Models 


In order to obtain a generic 6 DOF model, environment is also modeled. 
Gravity and atmosphere is modeled for simulations and derivations for wind effects 


are taken into consideration. 


For the gravity, WGS84 gravity model is used considering altitude of the 


vehicle above sea level [81]: 
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(4.28) 





JY) 
& =0.7803267714-[ 1+0.00193185138639- sin“ 4 | 


v1 —0.00669437999013-sin? 4 





g(h) = 8 | va] (4.29) 


where, 
A: Geographic latitude 


g,: Theoretical gravity 


R: Radius of the spherical earth 


h: Height of the vehicle above sea level 


For the atmosphere model, the 1962 International Standard Atmosphere or 
ISO 2553 is used given in [79]. Actually, since mid-course flight of a cruise missile 
at constant speed and altitude is considered, only wind model is required for the 


simulations. 


In 6 DOF simulations, winds and gusts alter the incidence angles and thus 
change the aerodynamic forces and moments. The incidence angles are calculated 
from the velocity vector of the vehicle’s center of mass with respect to air. To 


determine velocity of the vehicle with respect to air V,,,,,, wind vector Vp 1S 
subtracted from the geographic velocity v,,, considering the conventions given in 


Figure 103 [79]. 
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Vatair < YAE ~ Vairi E (4.30) 


—() Abe) (—(8) —(g) 

Valo =C us Ez) 

—( ) aa . 

Ve aa [-V,, “COS v, V “sin v, 0] (4.3 1) 
where, 


V,: Wind magnitude 


y,,: Wind direction from north 


In the simulations, a pre-specified wind profile with varying magnitude and 


direction is applied. Wind features will be given in the next section. 


4.1.2. Autopilots and Controls 


A variety of controllers are implemented for the mid-course flight of the 


cruise missile: 


1. Machhold control 


2. Roll position control 


3. Heading angle control with bank-to-turn autopilot 


4. Yaw stability augmentation 
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5. Altitude hold control with acceleration autopilot 


All autopilots are derived by classical pole placement techniques 


summarized in Zipfel [79]. 


4.1.2.1. Mach Hold Control 


Cruise missiles have to maintain Mach number under maneuvers and 
environmental effects. The thrust required F. to maintain a certain Mach number is 


equal to drag force projected onto the centerline of the turbine. Considering the 


turbine axis parallel to the body’s first axis, it is required that [79]: 


AO, 
cosa 


(4.32) 


Mach hold control loop is shown in Figure 105. The time constant 7, ofa 
generic turbojet engine is between 0.2 and 1.0 seconds. Gain G,, is calculated from 


second degree closed loop transfer function (TF) considering natural freguency and 


damping of the system as follows [79]: 


mV. 


mM 7 Mİ (4.33) 


where, 
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m: Mass of the missile 


V.: Sonic speed (V,=/yv-R-T ) 


¢: Damping of the closed loop TF 





Figure 105. Mach Hold Control Loop [79] 


4.1.2.2. Roll Position Control 


A dual feedback controller is built for roll position autopilot. The inner rate 
loop augments the aerodynamic damping and the outer loop executes the roll 
position command as shown in Figure 106. The transfer function between roll and 


aileron command is rather simple as follows [79]: 





p(s) to LL 5, 
da(s) s-LL, (134) 
where, 
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LL,,: Roll control derivative, LL ,,-(4Sb/1,)C, 


LL,: Roll damping derivative, LL, -(4Sb/7, )(b/2V ) C, 





Figure 106. Roll Rate and Position Feedback Loops [79] 


Using closed loop TF of the roll position against commanded roll position 


and setting parameters for second order TF, roll position autopilot gains K, and 


K, can be determined as follows [79]: 





(s) N K, al 
AA) s°+(K,-LL,,-LL,)s+K,-Ll,, SSS 
— 
26 onu, ro 
pen Bs (4.36) 
EE | 
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26 eyl, e LL p 


au (4.37) 


K 
di LL 





oa 


4.1.2.3. Heading Angle Control 


Heading changes of the cruise missiles are executed by roll control since 
bank-to-turn control is used. As the lift vector is banked, a horizontal force 
component generates a lateral acceleration that turns the velocity vector horizontally 


[79]: 


Heading angle tracker is built by wrapping a heading loop around the roll 
position autopilot as shown in Figure 107. Again, the pole placement technigue 


from root locus analysis is applied in order to determine heading gain K, using the 


open loop TF derived from the figure and eguation (4.35) [79]: 





NS 


Figure 107. Heading Angle Tracker Loop [79] 


2 


g Hs) gZ O, 
TF, (s)=K =>. =K >. = 4.38 
als) =K, Vs Bs) "V s*+ 26,30, 5 +O, s (138) 


Moll 
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y 
K, = O (1-67) (4.39) 


4.1.2.4. Yaw Stability Augmentation 


Dynamic stability of the vehicle in yaw plane is improved by the yaw rate 
damping loop. Yaw rate feedback loop is shown in Figure 108. The transfer 
function between yaw rate and rudder command is found from the linear 
perturbation eguations of the missile yaw plane which can be found in various 


references [79]: 





Figure 108. Yaw Rate Feedback Loop [79] 


res) _ ENa [s-¥p/V +(LN3/LNo İİ, G (5+2) 
Ör(s) —s?-(Y,/V+LN,)s+LN,+IN,Y, /V —s*+as+b 





(4.40) 


where, 


Y.: Dimensionalized derivatives of i for side force Y 
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LN,: Dimensionalized yawing moment of i 
Ys = (75/m)C;, > Ys, =(qS/m)C,, , 


IN, -(45b/1,)C,, , LN,, =(qSb/1)C,, , LN, =(qSb/1,)(b/2V )C,, l 


Using closed loop TF of the yaw rate against commanded yaw rate derived 
from Figure 108 and equation (4.40) and selecting the closed loop damping 


coefficient ¢., yaw rate gain K, is determined as follows [79]: 








r(s) - G, (s+z) N G,(s+z) 

Kis) s? +(a+K,G,)s+b+K,G,z  8°+26,,,0, s+0, (4.41) 
1 

K, = AĞ — 267.42) + Ya 2 eh -(@ 43,0) (4.42) 


o, =fb+K,Gz (4.43) 


4.1.2.5. Altitude Hold Control 


In order to build altitude hold autopilot, two feedback loops are wrapped 
around the normal acceleration autopilot with two gains G, and G, determining 
the dynamic response as shown in Figure 109. These gains are determined from the 
root locus analysis of the inner and outer altitude loops. Adaptive gain scheduling is 


not required for G, and G, , since altitude corridors are usually fixed and a constant 
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set of gains is sufficient. Actually, constant gain values give also good performance 


for terrain following (79). 














Figure 109. Altitude Hold Autopilot [79] 


4.1.2.6. Acceleration Autopilot 


In cruise missiles, the normal load factor plane generally contains an 
acceleration feedback loop. Guidance systems of cruise missiles like terrain 
following and obstacle avoidance require rapid response that only an acceleration 


autopilot can provide [79]. 


Acceleration autopilot loop is shown in Figure 110. The transfer function 
between pitch rate, normal acceleration and elevator command is found from the 
linear perturbation equations of the missile pitch plane which can be found in 


various references [79]: 
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Figure 110. Acceleration Autopilot Loop (79) 





M A 
qj | GE Ja M se | 
HE : PCI , Z (4.44) 


La: Dimensionalized derivative of lift force L with respect to a. 


M,: Dimensionalized pitching moment of i 
E; = (7S/m)C, , 


M, =(3$5c/1,)C, , Ms, =(3Sc/I,)C,, » M,= (3Sc/1, )(c/2V)C,,, , 


Using closed loop TF of the normal acceleration against commanded normal 


acceleration derived from Figure 110 and eguation (4.44) and selecting the closed 
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loop damping coefficient ¢ 


pilih? the natural frequency O pitch and a pole location p, 


acceleration autopilot gains are determined as follows [79]: 








ay(s) - M;VL,(G,s 4G, ) 
a(s) sV+s*(L,-MV+M,kV)+.. 
4.45 
-M L, -M ,V +M, kVL, +... SO 
KY + M 3,G,VL, 
M;kıL, +M ;,G,VL, 
O ap 
itch 
G, = = (4.46) 
a oe 
k li Oy, EDEM -že (4.47) 
2 M, pitch“ pitch q V : 
1 M,L M.L 
k > e 0 itc. +2 itc. © itch +M, + Z k, ee e) G, (4.48) 
1 LM, | pitch pitch“ p V V j; 


As it can be seen from the equations, the position feed-forward gain G, can 


not be determined from pole placement technique. Therefore, it must be determined 


from root locus analysis. Fixed value of gain G, determined from root locus gave 


sufficient performance for acceleration autopilot designed. 


4.1.3. Actuators 


An actuator is a device that actualizes steering inputs to motivators. These 


motivators are aileron, elevator and rudder for the cruise missile. For the 
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simulations, the response of the fin actuator is modeled by a second order TF as 
[79]: 





6i(s) O; 
= = 4.49 
ôi (s) s*+2C wet, S + O (543) 
where, 
Ol: Actual control surface deflection 


oi,: Fin command 
©, : Natural frequency of the actuator 


Ca: Damping ratio of the actuator 


Although the TF models only the linearized dynamics, fin deflection and fin 
rate limiters are included in the actuator model as two important nonlinearities 


which are shown in Figure 111. 








Figure 111. Second Order Actuator Model [79] 
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Actuator model derived is applied for all control channels (i.e. elevator, 
rudder and aileron) considering critical damping for the actuators by selecting large 
natural freguency values compared to autopilot natural freguencies. (For the 


simulations, &, =100 rad/s is selected for all channels.) 


4.1.4. Sensors 


Cruise missile model contain IMU and altimeter sensors. Accelerometers 
and gyros are used in the IMU of the missile. For height channel stability of the 
INS, barometric altimeter is used. For terrain clearance measurements in order to 


apply TAN algorithms, a radar altimeter is also reguired. 


Accelerometers are modeled considering strapdown INS architecture with 


random bias and noise, scale factor and misalignment as follows [79]: 





O FO 4 § FO (4.50) 
öf =V+(S,4+M,)- JO +W, (4.51) 
where, 


fc: Output of the accelerometers 


f: True acceleration values 


of: Accelerometer errors 


V, W: Random bias and white noise vector 
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S =| 0 asSF, 0 |: Scale factor error matrix 


M, = -aMA,, 0 aMA, |: Misalignment matrix 


Gyro errors are also modeled in the same manner [79]: 
Bi, = Pori +80 

600) ==+(5,+M,)-09+w, 

where, 

ON : Output of the gyros 

True angular rates 

060: Gyro errors 

E€, w,: Random drift and white noise vector 

S: Scale factor error matrix 


M,: Misalignment matrix 
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(4.52) 


(4.53) 


Barometric and radar altimeters are modeled with scale factor error and 
random noise. Actually, there exist detailed error models in the literature especially 
for barometric altimeters [82]. However, simple error models are sufficient for 
simulation purposes, since the effects of other errors can be neglected. Then error 


models for the altimeters become as follows [83]: 


ON is = Sharo Pa Viro (4.54) 
ON ise = S radan ` Padar * Vradar (4.55) 
where, 

öh: Altimeter errors 

S: Scale factors 

v: Measurement white noises 


4.2. Simulations 


In the previous section of the chapter, 6 DOF simulation models were 
developed. In this section, simulation results of a realistic cruise missile operation 


scenario with TAN will be discussed in detail. 


The operation scenario of the cruise missile model used in the simulations 


can be summarized as follows: 


1. 6 DOF cruise missile model with TAN simulations are performed 


with the loosely coupled architecture shown in Figure 101. 
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2. Simulation path is selected such that the missile had a constant 
heading through north for a constant time with large position errors. 
Then, TAN algorithms are applied during this time period in order to 
correct position errors. Next, the heading is corrected at discrete time 
intervals with the applied TAN algorithms in order to reach the 
desired waypoint. Simulation path used for the simulation is shown 


in Figure 112. Actually, selected terrain profile is a rough terrain. 





Na 
—İnitial Position - 
EN. a, 
ls Fi | jemi 40. 


sa PÄÄ 
Figure 112. Simulation Path Used for 6 DOF Cruise Missile Simulation 
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Large DTED grid size is used in acquisition mode for TAN 
algorithms and after first INS update, DTED grid size is decreased in 


tracking modes of operation. 


10 nmi/hr class INS sensor model parameters are used in order to 
model the INS used in the simulations. Initial errors are used for the 
INS in order to simulate cruise missile model, such that simulations 


began after about half an hour of operation. 


In order to accelerate simulation times, simulations are performed 
with constant aerodynamic parameters derived at the cruise Mach 
number. Hence, aerodynamic coefficients become function of angle 


of attacks only. 


Control strategy used for the simulations is as follows: 


a. Terrain following guidance is applied throughout the simulations. 
Altitude hold control is performed with pitch acceleration 


controller. 


b. Constant speed is achieved by Mach hold control loop throughout 


the operation. 


c. Heading angle tracker is used with bank-to-turn controller. 
Constant heading for initial period of operation is used. Then, 


heading updates are done at discrete times when INS is updated. 


d. Yaw stability is obtained with yaw stability augmentation 


controller. 
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Table 26. 6 DOF Cruise Missile Model Simulation Parameters 








Simulation Scenario 


1. 0—41 s. : Heading = 0° 

2. At4ls. : Initial heading correction 
for the desired waypoint 

3. At81s.: Second heading correction 





Controls Applied 


1. Terrain following with altitude hold 

2. Heading angle tracker with BTT 
control 

3. Yaw stability controller 

4. Machhold control 





Environmental Conditions 


1. ISO 2553 atmosphere model 
2. Pre-specified wind profile 
3. Universal gravity model 








Commanded Mach Number 0.7M 
Commanded Height 300 m AGL 
10 nmi/hr Class 


INS Ouality 


(Parameters are given in Table 5) 





Initial Horizontal INS Position Error 




















(Total Horizontal Position Error) oe 
Initial Horizontal INS Velocity Error 
Standard Deviation 0.5 m/s 
(Each axis, | o) 
Initial Vertical INS Position Error Si 
Standard Deviation (1 ©) 
Barometric Altimeter Standard 

ce 3m 
Deviation (1 ©) 
Radar Altimeter Standard Deviation 3 

m 

(10) 
Initial INS Attitude Errors Standard N 

ion 0.05 
Deviation (1 ©) 
TAN Algorithms Update Interval 0.5s 
INS Update Interval 40s 
TERCOM Update Interval 55s 
Initial DTED Grid Si 
TIA ae ee 11 x 11 (for PDAF and TSF) 
(Acquisition Mode) 











Tracking Mode DTED Grid Size 





5 x 5 (for PDAF and TSF) 
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Parameters used in the simulations are summarized in Table 26. Simulations 


are performed 


in real-time with these parameters for three INS position updates. 


Next, simulation results are presented. First, 6 DOF simulation results of the system 


are given as follows: 


Altitude vs. Time 
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Figure 113. Altitude vs. Time 
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Figure 114. Flight Profile over the Terrain 
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Horizontal Positions over Terrain Contours 
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Figure 115. Latitude vs. Longitude over Terrain Contours 
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Figure 116. Roll Rate vs. Time 
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Figure 117. Pitch and Yaw Rates vs. Time 
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Attitudes vs. Time 
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Figure 118. Attitudes vs. Time 
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Figure 119. Attitude Errors vs. Time 
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Figure 120. Total Velocity and Body Longitudinal Velocity vs. Time 
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Figure 121. Lateral and Vertical Body Velocities vs. Time 
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Figure 122. Angle of Attack and Side Slip Angle vs. Time 
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Figure 123. Wind Profile (Wind Velocity and Wind Heading) vs. Time 
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Figure 124. Mach vs. Time 
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Figure 125. Missile Mass vs. Time 
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Turbojet Thrust vs. Time 
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Figure 126. Turbojet Thrust vs. Time 
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Figure 127. Missile Heading and Commanded Heading vs. Time 
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Figure 128. Body Accelerations vs. Time 
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From the simulation results of the system, controls applied for mid-course 
guidance phase can be seen easily. In Figure 113 to Figure 115 flight profile of the 
missile is given in detail. Terrain following guidance and heading corrections for 
the reguired waypoint can be seen from the results. Commanded height is given as 
300 meters AGL and using acceleration controller, AGL height is kept between 300 
+ 150 meters which can be considered as an accepted range for the rough terrain 


considered. 


Attitude rates are shown in Figure 116 and Figure 117. Heading commands 
are given at 41 and 81 seconds of the simulation. Since BTT control is applied, roll 
and yaw rates change rapidly at these instants as expected. On the other hand, rapid 
changes in pitch rate are due to terrain following controller where normal 
acceleration controller is used in order to follow the rapid changes in height 
channel. Attitudes can be explained in the same manner which is shown in Figure 
118. Here, heading angle p changes due to missile bank angle ¢, since BTT 
controller is used. Moreover, pitch angle 0 changes between +10° due to terrain 
following. Finally attitude errors are presented in Figure 119 which are consistent 
with the INS error model results. Attitude and attitude rates in the controllers are 


used from INS error model. Hence, the system worked as in the actual applications. 


Missile total velocity and body velocity components are given in Figure 120 
and Figure 121. Total velocity is kept constant within a range due to Mach control 
loop. Changes in lateral body velocity are due to wind and rapid changes in vertical 
body velocity are due to terrain following control. Angle of attack & and side slip 
angle 4 are shown in Figure 122 and wind profile is given in Figure 123. Changes 
in 8 is due to wind profile where side wind components are considered for 
simulations. Finally Mach graph is presented in Figure 124 where Mach is 


controlled at 0.7 M within +0.002 M due to Mach control loop. 
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Changes in missile mass are given in Figure 125. Since Mach control is 
applied at constant altitude, constant mass flow rate, actually fuel flow rate is 
obtained as expected. Next, turbojet thrust is shown in Figure 126 where thrust 
changes between 750 — 2500 N in order to respond Mach control loop and terrain 
following controller. Then, missile heading versus commanded heading is given 
Figure 127. Here, missile response to heading commands can be seen well than the 
previous attitude plots. Heading corrections at 41 and 81 seconds of the simulation 
are applied considering the corrected missile position from TAN algorithms and 
waypoint position. It should be noted that INS updates from TAN algorithms are 
done at 40 seconds intervals. Hence, final position error decrease which will be 


discussed in the end of the chapter. 


Finally, body accelerations are presented in Figure 128 where longitudinal 
and lateral accelerations are small as expected due to Mach and BTT controls. 
Vertical acceleration component changes up to 3 g's due to acceleration controller 


in order to achieve terrain following. 


In the final section of the chapter, TAN results obtained from the 
simulations are presented. Here, original, updated INS and TAN algorithms? results 
are compared for both acguisition mode where initial position errors are huge for 
the first part of the simulation and for tracking mode. At the end of the chapter, final 


position errors are tabulated and discussed. 
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Figure 129. Northward Position Errors vs. Time (Acguisition Mode) 
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Figure 130. Eastward Position Errors vs. Time (Acquisition Mode) 
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Delta rN vs. Time (Tracking) 
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Figure 131. Northward Position Errors vs. Time (Tracking Mode) 
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Figure 132. Eastward Position Errors vs. Time (Tracking Mode) 
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Figure 133. Northward Velocity Errors vs. Time 
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Figure 134. Eastward Velocity Errors vs. Time 
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TAN results are compared for both acguisition modes and tracking modes of 
operation. In the first 40 seconds of the operation, TAN algorithms are applied for 
acguisition mode with a larger DTED grid size of 11x11 given in Table 26. Here, 
large initial position errors are corrected using PDAF and TSF algorithms. 
Moreover, the results are compared with TERCOM results. In Figure 129, 
northward position errors and in Figure 130, eastward position errors are shown. At 
40 seconds, PDAF and TSF corrections are done and INS is updated. On the other 
hand, TERCOM correction is done at 55 seconds. Both filters worked well in order 
to have errors less than 50 meters in both horizontal channels after a few seconds of 
operation. At the beginning of the TAN algorithms, past information is not used. 
Hence, filter results diverge for PDAF and TSF as it can be seen from the 


simulations. 


Then, tracking mode position errors are shown in Figure 131 and Figure 
132. Simulation results are given for 40 seconds to the end of the operation for 
tracking mode. In tracking mode, PDAF and TSF results are compared with SITAN 
results using a smaller DTED grid size of 5x5 given in Table 26. Here, INS is 
updated at 40 seconds intervals beginning from acquisition mode update at 40 
seconds. As it can be seen from the results, real-time PDAF results are better than 
other algorithms. TSF results follow INS error model as expected due to the 
algorithm properties. In TSF, index correction for the navigation system is done as 
in TERCOM. Hence, TSF errors are sometimes within the neighborhood of the 
actual navigation solution. However, errors are limited within the grid size of 
DTED Level 1 (i.e. -80 meters). Another interesting point is the divergence of TSF 
for a few seconds of operation at INS update times. This is due to lack of past 
measurement information. As older measurements come, TSF begins to follow INS 
errors. For PDAF, real-time corrections as in SITAN can be achieved since 


averages of the possible grid positions are considered throughout the operation. 


In Table 27, RMS errors and final position errors from the required 


waypoint is tabulated for tracking mode simulations. 
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Table 27. Tracking Mode Position Errors 























RMS Error [m] Final Position Error [m] 
PDAF 24 16 
TSF 60 75 
SITAN 36 24 
Updated INS 46 75 
Original INS 379 390 




















As it can be seen from Table 27, PDAF results are better than other 
algorithms, especially than SITAN. Hence, it can be concluded that real-time PDAF 
improves navigation performance. On the other hand, if TSF divergence is avoided 
for the initial few seconds of update intervals, similar results with the updated INS 
can be achieved. However, it is shown from the results that TSF works better for 
acquisition mode than tracking mode. For the simulations, loosely coupled 
architecture is used where INS is updated at discrete time intervals. Actually, 
simulation results show that tightly coupled architecture where INS is updated 
continuously can be used for real-time PDAF and SITAN in order to have better 


updated INS results. 


In the final part of the TAN simulations, horizontal position errors are 
presented in Figure 133 and Figure 134. TAN algorithms tend to correct northward 
velocity errors but not eastward velocity errors. Actually, with TAN algorithms 
velocity states can not be updated correctly, since measurements are only heights. 
However, if a detailed INS error model was used in TAN filters instead of simple 
error models derived in the previous chapter, better velocity errors could be 
obtained. In fact, TAN filter results can be improved by using better INS models in 
TAN models. 


259 


CHAPTER 5 


DISCUSSION AND CONCLUSION 


In the study, modern radar data association algorithms are implemented as 
new TAN algorithms which can be used with low-cost IMU’s. After performing a 
thorough survey of the literature on mid-course navigation of cruise missiles, study 
on modern radar data association algorithms and their implementations to TAN are 
done. Finally, performances of the designed navigation systems with the 
implemented TAN algorithms are examined in detail with the help of the 


simulations performed. 


In Chapter 1, theory about the study is given. Cruise missiles, cruise missile 
navigation performance and literature survey on TAN techniques are discussed. As 
it was discussed in detail in the chapter, the hearth of TAN is the algorithms used 
which fall into two general algorithmic categories of batch and recursive 
algorithms. Therefore, main research of the study is concentrated on TAN 


algorithms. 


From the literature survey, papers of Quintang, et al [40] and Dezert [43] 
gave inspiration for implementing modern data association algorithms to TAN in 
the Ph.D. study. Quintang, et al [40] propose a new TAN approach using PDAF to 
overcome irresolvable ambiguities in the correlation function used in TERCOM. 


The approach proposed is a batch algorithm which uses one of the modern radar 
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tracking algorithms. On the other hand, in the second paper, Dezert (43) proposes a 
new application of PDAF for improving the accuracy of autonomous strapdown 
INS. However, it is a real-time application of PDAF and relation with the former 
paper of Ouintang, et al [40] can be obtained where batch implementation of PDAF 
is used. Therefore, it was thought whether PDAF could be used as a TAN algorithm 


for real-time applications. 


TAN is a nonlinear estimation problem; since, terrain height information is 
used for navigation solution. Actually, TAN can be considered as a data association 
problem, especially for the acguisition operation mode where INS position errors 
are considerably large. From the literature survey as stated above, it has been 
thought that modern data association algorithms can be implemented for real-time 
TAN algorithms. Therefore, radar tracking, especially data association subject is 
investigated. At the end of Chapter 1, information about radar tracking technigues 


and possible implementations of radar data association algorithms to TAN is given. 


In Chapter 2, major TAN methods are investigated. First, INS errors of the 
cruise missiles and need for TAN systems are discussed. Then, major TAN methods 
including TERCOM, SITAN and VATAN are presented in detail. Fundamentals of 
the major methods are discussed in this chapter with simulations in order to make 


comparisons for the implemented TAN algorithms in the Ph.D. study. 


For the TERCOM process, several conclusions are achieved from the 


concept study and simulations performed. They are summarized as follows: 


1. Correlation algorithm is simple but not smart. Many calculations 
should be performed in order to have a position fix and navigation 


solutions can be obtained for rough and unigue terrains as expected. 
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It is thought that the algorithm was derived considering the 
capability of the computers of 1950's, performing only matrix 


calculations and simple mathematical operations. 


Physical meaning of MAD and MSD processes is the minimization 
of the area difference between the measured and the reference areas 
along the route of the missile. Actually, TERCOM process is 
actually a Maximum Likelihood Estimator (MLE) which uses “Least 
Sguares Estimation (LSE)” technigue. 


In the simulations, it was shown that MAD process shows better 
position fix than MSD process. For a terrain with small terrain height 
changes, MSD process neglects the small height difference terms and 
exaggerate the larger height difference terms. On the other hand, in 
MAD process absolute height difference terms are taken into account 


with same weights. 


The critical parameter for best terrain correlation is sigma-Z value of 
the area concerned where standard deviation of the point-to-point 
changes in terrain elevation (i.e. the slope) are calculated instead of 
sigma-T value where standard deviation of height of the area is 
calculated. In other words, the slopes of the area concerned are more 


critical than the roughness of the area for correlation. 


TERCOM process is independent of the target model where cruise 
missile is the target. Possible tracks for the missile are selected 
where tracks are the missile path formed by the terrain elevation file 
(DTED). Since, the target motion is not modeled; kinematical 


behavior of the system is not known. 
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For the SITAN process, conclusions achieved from the concept study and 


simulations performed is summarized as follows: 


1. SITAN is a recursive TAN technique which uses EKF unlike 
TERCOM which is a batch process. 


2. SITAN performance depends on the linearization of the terrain 
profiles since terrain slopes are reguired for the KF measurements. 
For large position errors, divergence can occur due to linearization 
errors in the EKF. In order to get rid of this, modified terrain 


linearization technigues and parallel KF structure are used. 


3. SITAN improves position errors for rough and mountainous terrain 
types. However, due to slope determination process in SITAN, 
solutions have sometimes serious jumps for mountainous terrain 
type. This can be explained by the severe slope changes in the 
mountainous terrain modeling. Therefore, linearization of the terrain 
profiles is very critical especially for mountainous terrains in 


SITAN. 


4. SITAN performance is better than both INS and terrain grids unlike 
TERCOM. In TERCOM, error can not be better than the terrain grid 


dimensions. 


5. SITAN performs better for smaller position errors due to terrain 
linearization. Due to this fact, for large initial position errors 


TERCOM or SITAN with parallel KF structure must be used. 
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In the last part of Chapter 2, VATAN is investigated. VATAN uses VA 
which is a maximum a posteriori (MAP) estimator that estimates a seguence of 
system states from a seguence of observation values. VA is actually a dynamic 
programming technigue for estimation which uses past information in data 
association problems. Therefore, it is thought that other data association algorithms 


can be used for TAN algorithms in the study. 


In Chapter 3, implementation of target tracking algorithms to TAN is 
presented. First, general information about modern target tracking algorithms are 
given. Next, PDAF and TSF data association algorithms and their general 
implementations are investigated. Then, PDAF and TSF implementations to TAN 
are presented. At the end of the chapter, a simple simulation model is developed for 
the mid-course flight of the cruise missile. Finally, simulations are performed with 
the implemented TAN algorithms and the results are compared with the major TAN 


methods. 


The advantages of the PDA and TS approach implemented for TAN solution 


can be summarized as follows: 


1. Real-time TAN solution can be obtained with a single PDA filter or 
parallel TS filters. 


2. PDA and TS filters can be used for both batch and recursive TAN 
solution. For batch solution, larger grid size is selected for navigation 
solution. For recursive solution, horizontal positions are calculated 


recursively in relatively small DTED grids. 


3. Since past measurements are taken into account, smoothing of the 


measurements in the filter is achieved which decreases errors. 
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4. Since INS error model is used for navigation solution, application of 


the filters is simple and the filters are linear. 


Batch size of the DTED area concerned can be changed. Both larger 
DTED areas for acguisition mode or smaller DTED areas for 


tracking modes can be selected using the same PDA filter. 


TSF gives solutions for various tracks selected. Actually, all tracks 
converge to the same index of the DTED grid (i.e. solution grid). 
However, for smooth terrains, there exist more than one position 
solution index and the tracks can be investigated separately in order 


to give more than one but finite number of navigation solutions. 


TSF approach is original when compared with other papers. However, PDA 


approach for TAN is found in the literature. The difference of the PDA algorithm 


developed from Oingtang, et al [40] is summarized as follows: 


1. 


In the paper of Oingtang, et al [40], TAN using PDAF was 
investigated for the batch algorithm. The motion of the vehicle is not 


modeled. 


In the paper of Oingtang, et al [40], performance of the TAN using 
PDA and TERCOM has been compared. It is stated that PDA was 
used in order to improve the performance of TAN compared to 


TERCOM. 


In the Ph.D. study, real-time PDAF implementation is done. By 


using the error model of the INS used in the vehicle, system 
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dynamics is modeled. Using PDAF, error states of the system are 


estimated recursively. 


4. In the Ph.D. study, PDAF equations are directly implemented for the 
TAN solution. Association probabilities obtained from height 
difference measurements for each element of the DTED grid 
concerned are used for position updates, considering the index of the 


DTED grid. 


In the last part of Chapter 3, simulations are performed for both acguisition 
and tracking modes of operation considering a small period of time (i.e. 100 
seconds for tracking mode) operation in order to visualize the performance of the 
implemented TAN algorithms. Simulations are performed for rough, smooth and 
mountainous terrain types. Moreover, effects of using different DTED types and 


DTED grid sizes are also investigated. 


In tracking mode, it is seen that better results than SITAN are obtained for 
rough and mountainous terrain types. TSF and PDAF results are considerably good, 
since 1f navigation solution does not exist, the filters follow INS error model which 
is actually a desired feature. From the Monte Carlo simulations, position RMS 
errors of the TSF and PDAF algorithms become less than 50 meters for 


mountainous terrains; in other words, a decreased navigation error is obtained. 


In acguisition mode, it is seen that similar results with TERCOM are 
obtained for rough terrain type. Critical point in the acguisition mode simulations is 
the percentage of false fix for acguisition mode. In TERCOM, since INS error 
model is not considered for correlation process, there is always a probability of false 


fix in the position solutions. In the Monte Carlo simulations performed, a few false 
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position fixes occurred for TERCOM for rough terrain type. On the other hand, 
with PDAF and TSF no false position fixes occurred. 


Simulation results also show that TAN algorithms work with DTED Level 2 
for acguisition and tracking modes. However, DTED Level 1 results seem to be 
better than DTED Level 2 results. Unfortunately, there were not sufficient DTED 
Level 2 maps for simulations in order to compare simulation results in detail. 
Actually, vehicle velocity directly influences TAN performance. For cruise 
missiles, DTED Level 1 maps are sufficient for mid-course flight navigation 
solution where INS position fixes less than 50 meters can be obtained. For faster 
vehicles like cruise missiles, rapid changes in the terrain profile as in DTED Level 2 
decreases TAN performance. As a result of this, use of DTED Level 1 maps for 


TAN acguisition mode seems to perform better solutions. 


At the end of Chapter 3, simulations are done with various DTED grid sizes 
for PDAF and TSF. Here same initial position errors are taken for simulations in 
tracking mode along rough terrain. It is seen that changing grid size for the same 
initial position errors for TSF slightly changes simulation results. Again errors are 
bounded and limited with the related grid solution. However, selecting larger grid 
sizes for PDAF solutions generally increase position errors. This is due to PDA 
procedure where weighted averages of the all grid points are taken into account for 
navigation solution. Therefore, it can be concluded that PDAF DTED grid size 


should be selected in accordance with the position errors. 


In Chapter 4, case studies are performed for a cruise missile model with the 
help of the 6 DOF simulation tool developed. The simulation tool developed for the 
cruise missile is capable of performing full mid-course flight simulation of the 
cruise missile modeled. Actually, a generic simulation tool applicable to all air 


vehicles is considered except for guidance methods applied. 
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A variety of controllers are implemented for the mid-course flight of the 
cruise missile; Mach hold control, roll position control, heading angle control with 
bank-to-turn autopilot, yaw stability augmentation and altitude hold control with 
acceleration autopilot. All autopilots are derived by classical pole placement 


technigues summarized in Zipfel [79]. 


Then, simulations are performed with PDAF and TSF TAN models with 
actual flight conditions. Finally, simulation results are compared with major TAN 
algorithms considering other flight parameters of the cruise missile model. From the 
simulation results of the system, controls applied for mid-course guidance phase are 


clearly observed. 


From the simulations performed in Chapter 4, better results are obtained for 
PDAF than other algorithms, especially than SITAN. Hence, it can be concluded 
that real-time PDAF improves navigation performance. However, simulation results 


also show that TSF works better for acguisition mode than tracking mode. 


Several conclusions are achieved from the implemented PDAF and TSF 
algorithms from the simulations performed in Chapter 3 and Chapter 4. The 


advantages of the new algorithms proposed can be summarized as follows: 


1. Real-time TAN solution can be obtained with a single PDA filter. 
Since past measurements are taken into account, by changing the 
buffer size of the measurements the filter, measurements are 


smoothed. 


2. Real-time TAN solution can be obtained with a single TSF structure. 
However, TSF operations are more complex than SITAN. On the 
other hand, in TSF, more than one track is selected in order to 
determine navigation solution. Hence, probability of false fix 


decreases unlike in TERCOM. 
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Real-time TAN solution is obtained by considering horizontal 
position errors of DTED used in real-time PDA filter and TSF. 
Hence, horizontal position states are added to the Kalman filters used 


in PDAF and TSF. 


. Application of the filters is simple and the filters are linear, since 


INS error model is used. 


Batch size of the DTED area concerned can be changed independent 
of the model used. Both larger DTED areas for acguisition mode or 
smaller DTED areas for tracking modes can be selected using the 


same filters. 


Results of the filters are good for both recursive and batch 
algorithms. For tracking mode, position RMS error is less than 50 
meters for DTED Level 1. Moreover, PDAF shows stable response. 
For smooth terrains where no navigation solution exists, PDAF 


follows the INS error model which is actually a desired feature. 


TSF can be considered as a real-time TERCOM process for large 
position errors, i.e. large DTED batch size. Possibility of false 
position fixes decrease with TSF when compared with TERCOM. 
On the other hand, for small position errors, decreasing the 
weighting factor of the past measurements for TSF, better real-time 
solutions can be obtained. However, real-time results of TSF follows 


INS error model unlike PDAF where position errors decrease much. 
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Contributions of the Study: 


In order to declare the contributions of the study, disadvantages of the 


present TAN algorithms should be overviewed again. Several disadvantages of the 


TAN algorithms can be summarized as follows: 


TAN reguires terrain information for real-time navigation solution 
and the dynamics of the system is highly nonlinear which need 


considerable calculation work. 


Real-time application for the TAN solution is generally impractical 
for high velocity vehicles like cruise missiles due nonlinear 


characteristics of the system. 


In SITAN, terrain linearization and terrain slopes are reguired in 
order to apply extended Kalman filter eguations which are actually 


critical stages for TAN solution. 


TERCOM is a batch process and it is independent of the target 
model where cruise missile is the target. Since, the target motion is 
not modeled; kinematical behavior of the system is not known and 
possibility of false position fixes increase especially for terrains with 


similar height profiles. 


Then, the contributions of the study can be summarized as follows: 


1. 


Modern radar data association algorithms are implemented as new 


TAN algorithms which can be used with low-cost IMU’s. 
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. Acquisition mode performance of the TAN algorithms is improved 
when compared with TERCOM. Probability of false fix decreases 
with the implementation of PDAF and TSF for TAN. 


Tracking mode performance of the TAN algorithms is improved 
when compared with SITAN especially with the implementation of 
PDAF. 


. Application of the filters is simple and the filters are linear, since 


INS error model is used for position updates. 


. No linearization for terrain is reguired for the implemented 
algorithms. DTED files can be used directly without any prior work 


for operation. 


. Implemented algorithms can be applied to existing systems with the 


use of the new micro-processors with relatively low costs. 


Future Work: 


As the future work, implemented TAN algorithms can be used for tightly 


coupled integration architecture. Actually, simulation results show that tightly 


coupled architecture where INS is updated continuously can be used for real-time 


PDAF and SITAN in order to have better updated INS results. 


Another future work can be the improvement of the implemented filters. 


Filter system and measurement models are constructed considering simple INS 


models with position and velocity errors only. If a detailed INS error model was 


used in the implemented filters, better results could be obtained than the simulation 
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results in Chapter 4. Finally, a hybrid filter algorithm can be implemented which 
uses both PDAF and TSF in order to have superior results than stand alone PDAF 
and TSF algorithms. 


272 


REFERENCES 


[1]. Wikipedia, the free Encyclopedia, “Cruise Missile”, 
http://www.wikipedia.org/w/wiki.phtml?search=cruise+missile&g0=Go, 


Last accessed in May 2003. 


[2]. Siouris, G., M., Missile Guidance and Control Systems, Springer-Verlag 


New York, Inc., pp. 521-522, 2004. 


[3]. Metzger, J., Wendel, J., Trommer, G., F., “Hybrid Terrain Referenced 
Navigation System using a Bank of Kalman Filters and a Comparison 
Technigue”, AlAA, Guidance, Navigation, and Control Conference and 


Exhibit, Providence, Rhode Island, August 2004. 


[4]. Federation of American Scientists, "Cruise Missiles”, 


http://www.fas.org/nuke/intro/cm, Last accessed in May 2003. 


[5]. Federation of American Scientists, “BGM-109 Tomahawk”, 
http://www.fas.org/man/dod-101/sys/smart/bgm-109.htm, Last accessed in 


May 2003. 


273 


[6]. 


[7]. 


[8]. 


[9]. 


[10]. 


[11]. 


Federation of American Scientists, “4GM-86 Air-Launched Cruise Missile 
[ALCM] ”, http://www.fas.org/nuke/guide/usa/bomber/ alcm.htm, Last 


accessed in May 2003. 


CDISS, the Centre for Defence and International Security Studies, “Cruise 
Missiles; Key Technologies: An Overview”, 


http://www.cdiss.org/tabtechs.htm, Last accessed in May 2003. 


DARPA Special Projects Office, “Low Cost Cruise Missile Defense 
(LCCMD) ”, http://www.darpa.mil/spo/Programs/lowcostcruisemissile- 


defense.htm, Last accessed in May 2003. 


DARPA Special Projects Office, “LCCMD Side Show”, 
http://www.darpa.mil/spo/Programs/LCCMD slides/ 


LCCMD slide show.pdf, Last accessed in May 2003. 


Aardvark, “4 DIY Cruise Missile”, http://www.interestingprojects.com/ 


cruisemissile/, Last accessed in May 2003. 


Johnson, N., Tang, W., Howell, G., “Terrain Aided Navigation Using 
Maximum A Posteriori Estimation”, IEEE, Position Location and 


Navigation Symposium, 1990. 


274 


[12]. 


[13]. 


[14]. 


[15]. 


[16]. 


[17]. 


[18]. 


Hostetler, L., D., Andreas, R., D., “Nonlinear Kalman Filtering Techniques 
Jor Terrain-Aided Navigation”, IEEE Transactions on Automatic Control, 


Vol. AC-28, No. 3, March 1983. 


Boozer, D., D., “Terrain Referenced Navigation”, AGARD-AG-331, 


Aerospace Navigation Systems, pp. 152-157, June 1995. 


The MathWorks Inc., “MATLAB Help, The Language of Technical 


Computing”, Version 7.1.0.246 (R14SP3), 1984-2005. 


MIL-PRF-89020A, “Performance Specification Digital Terrain Elevation 


Data (DTED)”, US Department of Defense, 1996. 


“Sayısal Harita ve Coğrafi Bilgi Sistemi Kurs Notları”, T.C. M.S.B. Harita 
Genel Komutanlığı, Tasnif Dışı, TÜBİTAK-SAGE Kütüphanesi, CPP733, 


2003, in Turkish. 


Henley, A., J., “Terrain Aided Navigation — Current Status, Technigues for 
Flat Terrain and Reference Data Reguirements”, IEEE, Position Location 


and Navigation Symposium, 1990. 


Nielson, J., T., “CALCM — The Untold Story of the Weapon Used to Start 
the Gulf War”, IEEE Aerospace and Electronics Systems Magazine, July 


1994. 


275 


[19]. 


[20]. 


[21]. 


[22]. 


[23]. 


[24]. 


[25]. 


Hicks, S., "Advance Cruise Missile Guidance System Description”, 


Proceedings of the IEEE Aerospace and Electronics Conference, 1993. 


Bennett, P., J., “The Use of Digital Map Data for Airborne Operations”, 


TEE Colloguium on Serious Low Flying, Ref. No. 1998/223. 


Chen, Z., Yu, P., “Model Study for Terrain Aided Navigation Systems”, 
Proceedings of the IEEE International Symposium on Industrial Electronics, 


1992. 


Wang, W., Chen, Z., “Error Model Identification on Digital Map of TAN 
System Based on EKF” , Proceedings of the IEEE International Conference 


on Industrial Technology, 1994. 


Yu, P., Chen, Z., Hung, J., C., “Performance Evaluation of Six Terrain 
Stochastic Linearization Techniques for TAN”, Proceedings of the IEEE 


National Aerospace and Electronics Conference, 1991. 


Paris, S., Le Cadre, J., “Planification for Terrain-Aided Navigation”, IEEE, 
Proceedings of the Fifth International Conference on Information Fusion, 


2002. 


McFarland, M., B., Zachery, R., A., Taylor, B., K., “Motion Planning for 
Reduced Observability of Autonomous Aerial Vehicles”, Proceedings of the 


IEEE International Conference on Control Applications, 1999. 


276 


[26]. 


[27]. 


[28]. 


[29]. 


[30]. 


[31]. 


[32]. 


Bar-Gill, A., Ben-Ezra, P., Bar-Itzhack, I., Y., “Improvement of Terrain- 
Aided Navigation via Trajectory Optimization”, IEEE, Transactions on 


Control Systems Technology, Vol. 2, No. 4, December 1994. 


Li, D., Zhou, D., Hu, Z., Hu, H., “Optimal Preview Control Applied to 
Terrain Following Flight”, Proceedings of the 40th IEEE Conference on 


Decision and Control, 2001. 


Baird, C., A., Snyder, F., B., “Terrain-Aided Altitude Computations on the 


AFTI/F-16”, IEEE, Position Location and Navigation Symposium, 1990. 


Hollowell, J., “Heli/SITAN: A Terrain Referenced Navigation Algorithm for 


Helicopters”, IEEE, Position Location and Navigation Symposium, 1990. 


Nordlund, P., Gustafsson, F., “Recursive Estimation of Three-Dimensional 
Aircraft Position Using Terrain-Aided Positioning”, IEEE International 


Conference on Acoustics, Speech, and Signal Processing, 2002. 


Bergman, N., Ljung, L., Gustafsson, F., “Point-Mass Filter and Cramer- 
Rao Bound for Terrain-Aided Navigation”, Proceedings of the 36th IEEE 


Conference on Decision and Control, 1997. 


Newman, P., Durrant-Whyte, H., “Using Sonar in Terrain-Aided 
Underwater Navigation”, IEEE International Conference on Robotics and 


Automation, 1998. 


277 


(331. 


[34]. 


[35]. 


[36]. 


[37]. 


[38]. 


[39]. 


Madhavan, R., Durrant-Whyte, H., Dissanayake, G., "Natural Landmark- 
Based Autonomous Navigation Using Curvature Scale Space”, IEEE 


International Conference on Robotics and Automation, 2002. 


Bruder, S., B., H., Wedeward, K., “Terrain Aided INS Robot Navigation: A 
Deferred Decision Making Approach”, IEEE, 42nd Midwest Symposium on 


Circuits and Systems, 1999. 


Morisue, F., Ikeda, K., “Evaluation of Map-Matching Technigues”, IEEE 


Vehicle Navigation and Information Systems Conference, 1989. 


McLellan, J., F., Schleppe, J., “GPS/Barometry Height-Aided Positioning 


System”, IEEE, Position Location and Navigation Symposium, 1994. 


Siouris, G., M., Missile Guidance and Control Systems, Springer-Verlag 


New York, Inc., pp. 551-576, 2004. 


Erhui, W., Guohua, G., Lincheng, S., Wensn, C., “A Probability-Based 
Terrain-Aided Navigation Approach and Its Relative Terrain Navigability 
Analysis”, Proceedings of the IEEE International Conference on Industrial 


Technology, 1996. 


Zhou, H., Zhang, C., “Terrain Aided Navigation Based on Computer 
Vision”, Proceedings of the 4th World Congress on Intelligent Control and 


Automation, 2002. 


278 


[40]. 


[41]. 


[42]. 


[43]. 


[44]. 


[45]. 


Oingtang, F., Lincheng, S., Wenseng, C., "Terrain Aided Navigation Using 
PDAF”, Proceedings of the 2003 IEEE, International Conference on 


Robotics, Intelligent Systems and Signal Processing, China, 2003. 


Pei, Y., Chen, Z., Hung, J., C., "BITAN-II: An Improved Terrain Aided 
Navigation Algorithm”, Proceedings of the 1996 IEEE IECON 22 Industrial 


Electronics, Control, and Instrumentation, 1996. 


Enns, R., Morrell, D., “Terrain-Aided Navigation Using the Viterbi 
Algorithm”, AIAA, Journal of Guidance, Control, and Dynamics, Vol. 18, 


No. 6, November-December 1995. 


Dezert, J., “Improvement of Strapdown Inertial Navigation using PDAF”, 
IEEE Transactions on Aerospace and Electronic Systems, Vol. 35, No. 3, 


July 1999. 


Hongbo, X., Yan, T., JianZhong, S., Jinwen, T., Jian, L., “Terrain Matching 
Based on Imaging Laser Radar”, IEEE 6th International Conference on 


Signal Processing, 2002. 


Bevington, J., E., Marttila, C., A., “Precision Aided Inertial Navigation 
Using SAR and Digital Map Data”, IEEE, Position Location and Navigation 


Symposium, 1990. 


279 


[46]. 


[47]. 


[48]. 


[49]. 


[50]. 


[51]. 


[52]. 


Chan, L., C., Snyder, F., B., “System for Correlation and Recognition of 
Terrain Elevation”, United States Patent, Patent No: 4,584,646, Date of 


Patent: Apr. 22, 1986. 


Baird, C., A., “Map-Aided Navigation System Employing TERCOM-SITAN 
Signal Processing”, United States Patent, Patent No: 4,829,304, Date of 


Patent: May 9, 1989. 


Lerche, H., “Navigation of Aircraft by Correlation ”, United States Patent, 


Patent No: 4,910,674, Date of Patent: Mar. 20, 1990. 


Raymer, K., A., McGuffin, J., T., “Terrain Referenced Navigation — Schuler 
Cycle Error Reduction Method and Apparatus”, United States Patent, Patent 


No: 5,450,345, Date of Patent: Sep. 12, 1995. 


Goebel, R., H., Fogle, D., A., Torretta, D., C., Panagos, P., Heffern, P., A., 
“Terrain Correlation System”, United States Patent, Patent No: 6,218,980, 


Date of Patent: Apr. 17, 2001. 


Uhlmann, J., “Introduction to the Algorithmics of Data Association in 
Multiple-Target Tracking”, Handbook of Multisensor Data Fusion, CRC 


Press LLC, 2001. 


Kirubarajan, T., Bar-Shalom, Y., “Target Tracking Using Probabilistic 


Data Association Based Techniques with Applications to Sonar, Radar, and 


280 


EO Sensors”, Handbook of Multisensor Data Fusion, CRC Press LLC, 


2001. 


[53]. Stieler, B., “Inertial Navigation”, AGARD-AG-331, Aerospace Navigation 


Systems, p. 45, June 1995. 


[54] Rogers, R., M., Applied Mathematics in Integrated Navigation Systems, 
AIAA Education Series, ISBN: 1-56347-397-6, Reston, VA, pp. 273-282, 


2000. 


[55] Bar-Itzhack,I., Y., Berman, N., “Control Theoretic Approach to Inertial 
Navigation Systems”, ALAA Journal of Guidance, Vol. 13, pp. 237-245, 


May-June 1988. 


[56] Scherzinger, B., M., Reid, D., B., “Modified Strapdown Inertial Navigator 
Error Models”, Proceedings of the 1994 IEEE PLANS, Las Vegas NY, 


April 1994. 


[57]. Phillips, R., E., Schmidt, G., T., “GPS/INS Integration”, AGARD-LS-207, 
System Implications and Innovative Applications of Satellite Navigation, p. 


9-12, June 1996. 


[58]. The MathWorks Inc., “Simulink, Dynamic System Simulation for 


MATLAB”, Version 6.3 (R14SP3), 1990-2005. 


281 


[59] 


[60]. 


[61]. 


[62]. 


[63] 


[64]. 


[65]. 


[66]. 


Kayton, M., Fried, W. R., ed., Avionics Navigation Systems, John Wiley & 


Sons, pp. 317-318, 1969. 


The MathWorks Inc., “MATLAB, The Language of Technical Computing”, 


Version 7.1.0.246 (R14SP3), 1984-2005. 


Newman, D., “OziExplorer, GPS Mapping Software”, Version 3.95.3.g2, 


D&L Software Pty Ltd., Australia, 2004. 


Newman, D., “OziExplorer3D”, Version 1.07b, D&L Software Pty Ltd., 


Australia, 2004. 


Gelb, A., ed., Applied Optimal Estimation, The M.I.T. Press, pp. 182-189, 


Eleventh Printing, 1989. 


Chong, C., Garren, D., Grayson, T., P., “Ground Target Tracking — A 


Historical Perspective”, IEEE, Aerospace Conference Proceedings, 2000. 


Durrant-Whyte, H., “Introduction to Estimation and Data Fusion”, Lecture 
Notes, Australian Centre for Field Robotics, The University of Sydney, 


2003. 


Chang, K., “Syst680: Principles of C31 #4”, Lecture Notes, Dept. of 


Systems Engineering, George Mason University, 2004. 


282 


(671. 


[68]. 


[69] 


[70]. 


[71]. 


[72]. 


[73]. 


Maksarov, D., Durrant-Whyte, H., “Mobile Vehicle Navigation in Unknown 
Environments: A Multiple Hypothesis Approach”, IEE Proc.-Control Theory 


Appl., Vol. 142, No. 4, July 1995. 


Stone, L., D., “4 Bayesian Approach to Multiple-Target Tracking”, 


Handbook of Multisensor Data Fusion, CRC Press LLC, 2001. 


Blackman, S., Popoli, R., Design and Analysis of Modern Tracking Systems, 


Artech House, ISBN No: 1-58053-006-0, pp. 360-402, 1999. 


Reid, D., B., “An Algorithm for Tracking Multiple Targets”, IEEE 


Transactions on Automatic Control, Vol. AC-24, No. 6, Dec. 1979. 


Kurien, T., “Issues in the Design of Practical Multitarget Tracking 
Algorithms”, Multitarget-Multisensor Tracking: Advanced Applications, Y. 


Bar-Shalom (Ed.), Norwood, MA: Artech House, 1990, Chap. 3. 


Smith, P., B., Buechler, G., “4 Branching Algorithm for Discriminating and 
Tracking Multiple Objects”, IEEE Transactions on Automatic Control, pp. 


101-104, Feb. 1975. 


Stieler, B., “Inertial Navigation”, AGARD-AG-331, Aerospace Navigation 


Systems, p. 100, June 1995. 


283 


[74]. 


[75]. 


[76]. 


[77]. 


[78]. 


[79] 


[80]. 


[81]. 


Kumar, M., “Coordinate Frames”, AGARD-AG-331, Aerospace 


Navigation Systems, pp. 7-41, June 1995. 


The MathWorks Inc., “MATLAB, Mapping Toolbox”, Version 2.2 


(R14SP3), 1984-2005. 


The MathWorks Inc., “Simulink, Signal Processing Blockset”, Version 6.2 


(R14SP3), 1990-2005. 


“Microdem/ TerraBase II Software”, PETMAR Trilobite Breeding Ranch, 


Version 14.45.15.2, March 2003. 


Özgören, M. K., “ME 502 Advanced Dynamics Lecture Notes”, Middle East 


Technical University, 1997. (Unpublished) 


Zipfel, P., H., Modeling and Simulation of Aerospace Vehicle Dynamics, 


AIAA Education Series, ISBN: 1-56347-456-5, Reston, VA, 2000. 


McDonnell Douglas Astronautics Company, “The USAF Stability and 
Control Digital DATCOM, Volume I, Users Manual”, Technical Report, 
AFFDL-TR-79-3032, Vol. I, St. Louis Division, St. Louis, Missouri, April 


1979. 


“International Gravity Formula”, http://geophysics.ou.edu/solid earth/ 


notes/potential/igf.htm, Last accessed in March 2006. 


284 


(821. 


(831. 


Lewantowicz, Z., H., Paschall, R., N., “Deep Integration of GPS, INS, SAR, 
and Other Sensor Information”, AGARD-AG-331, 4erospace Navigation 


Systems, p. 239, June 1995. 


Liang, D., F., “An Overview of a Generic Multi-Sensor Integrated 
Navigation System Design”, AGARD-AG-331, Aerospace Navigation 


Systems, pp. 219-220, June 1995. 


285 


APPENDIX 
PROBABILISTIC DATA ASSOCIATION EQUATIONS [52] 
The PDA algorithm calculates in real-time the probability that each 


validated measurement is attributable to the target of interest. This probabilistic 


(Bayesian) information is used in a tracking filter, the PDA filter (PDAF), which 


accounts for the measurement origin uncertainty. 


Past Measurement Information: 


The PDAF uses a decomposition of the estimation with respect to the origin 


of each element of the latest set of validated measurements, denoted as: 


ZK) = fz (A.1) 
where; 


z,(k): i’th validated measurement, 


m(k): Number of measurements in the validation region at time k. 


The cumulative set (sequence) of measurements is: 
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Z*-(2(0)., (A2) 


Measurement Validation: 


From the Gaussian assumption, the validation region is the elliptical region: 


V(k,7)=|Z:[2-2k]k-D] -S7 [2 -20k kD] <7} (A.3) 


where; 


y: Gate threshold, 
S(k) = H(k)-P(k|k-1)-H(k)' + R(k) (A.4) 


S(k): Covariance of the innovation corresponding to the true measurement. 


The volume of the validation region given in eguation (A.3) is: 


V(k)=c, - jä 


y Sİ =c, 7" İSA (A.5) 





where the coefficient c, depends on the dimension of the measurement (it 
is the volume of the n -dimensional unit hyper sphere: c, =2, c,=7, c,;=4n/3, 


etc.) 
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The State Estimation: 


In view of the assumptions listed, the association events; 


96) = [ is the target originated measurement} i=1,...,m(k) 


{none of the measurements is target originated} i=0 (A.6) 


are mutually exclusive and exhaustive for m(k)>1. 


Using the total probability theorem with regard to the above events, the 


conditional mean of the state at time k can be written as; 
R(k | k) = E| x(k)| Z* | 


m(k) 


He) = 3 E| x(k) |0,00),2" Paz) 


m(k) 


X(k| k)= 2421) (k) (A.7) 


where, x,(k|k) is the updated state conditioned on the event that the i’th 


validated measurement is correct, and; 
B,(k) = P{O(k)|Z"| (A.8) 


is the conditional probability of this event; the association probability, 


obtained from the PDA procedure presented in the next subsection. 


The estimate conditioned on measurement 7 being correct is; 
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EKİ) (k |k-D)+K(k)-v(k) i=1,...,m(k) (A.9) 


where the corresponding innovation is; 


v.(k) = z (k) -2(k|k-1) (A.10) 


The gain K(k) is the same as in the standard Kalman filter; 


K(k) = P(k|k-1)-H(k)’ -S(k)*' (A.11) 


since, conditioned on 0,(k) , there is no measurement origin uncertainty. 


For i=0 (i.e. if none of the measurements is correct) or m(k) =0 (i.e. there 


is no validated measurement); 


£,(k | k) = X(k|k-1) (A.12) 


The State and Covariance Update: 


Combining equations (A.9) and (A.12) into equation (A.7) yields the state 
update equation of the PDAF; 


X(k | k) = x(k|k-1)+ K(k)-v(k) (A.13) 


where the combined innovation is; 
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m(k) 


vk) = 2) ARAC) (A.14) 


The covariance associated with the updated state is; 

P(k | k) = B,(k)-P(k|k-1)+[1—8,(6)]-P*(Kk|k)+ PCR) (A.15) 
where the covariance of the state updated with the correct measurement is; 
P*(k | k) = P(k|k-1)—K(k)-S(k)- K(k)’ (A.16) 


and the spread of the innovations term (similar to the spread of the means 


term in a mixture) is; 


b mk) 
P(k)=K(k)- | 2 Bk) -vi (k) -v(k)' — v(k) ON K(k) (A.17) 


The Prediction Eguations: 


The prediction of the state and measurement to k+1 is done as in the 


standard filter, i.e.; 
X(k+1|k)=DP(k)-X(k | k) (A.18) 


2(k+1|k)=H(k+1)-2(k+1|k) (A.19) 
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The covariance of the predicted state is, similarly; 
P(k+1|k)=@(k)- P(k| k)-®(k)’ + O(k) (A.20) 
where P(k |k) is given by equation (A.15). 


The innovation covariance (for the correct measurement) is, again, as in the 


standard filter; 


S(k+1)= H(k+1)-P(k+1|k)-H(k+1)" + R(k+1) (A.21) 


The Probabilistic Data Association: 


To evaluate the association probabilities, the conditioning is broken down 
into the past data Z* and the latest data Z(k). A probabilistic inference can be 


made on both the number of measurements in the validation region (from the clutter 


density, if known) and on their location, expressed as; 
B,(k) = PG] Z"} PİİZ, mk), 2°" (A.22) 


Using Bayes’ formula, the above is rewritten as; 


f(t) =~ p[ Zam.) p {00 mZ) 
i=0,...,m(k) 


(A.23) 
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The joint density of the validated measurements conditioned on 4 (k),i #0, 


is the product of; 


e The (assumed) Gaussian PDF of the correct (target-originated) 


measurements; 


* The PDF of the incorrect measurements, which are assumed to be uniform 


in the validation region whose volume V(k) is given in eguation (A.5). 


The PDF of the correct measurement (with the P, factor that accounts for 


restricting the normal density to the validation gate) is; 


p| 2,00] 0,06), m(k), 2 | = pg!-N [2,06] | 2(k lk —-1),5(6)] 


= po -N[v,(k) |0,5(6)] 


(4.24) 


The PDF from eguation (A.23) is then; 


p| Z06)] 8.00, mk), Zİ |EN. (4.25) 


The probabilities of the association events conditioned only on the number 


of validated measurements are; 


V (km. P, Ny) | o, S(k)] i=l,...,m(k) 


A.26 
Vio i=0 ( ) 


pİZ O(k),m(k), Z |= | 
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where 4,(m) is the probability mass function (PMF) of the number of false 


measurements (false alarms or clutter) in the validation region. 


Two models can be used for the PMF £,(m) in a volume of interest V : 
1. A Poisson model with a certain spatial density 4; 


mme E (4.27) 
m! 


2. A diffuse prior model; 


Hp (m) = H;(m-1)=ö (A.28) 
where the constant 6 is irrelevant since it cancels out. 
Using the (parametric) Poisson model in equation (A.26) yields; 
P,-P,-[P,-P,-m(k)+(1-P,-P,)-4-V(k)]' i=1,....m(k) 
lm T | i (A.29) 
(1-P,-P,)-2-Vik)-[P, Em) AP, EŞ AVİ i=0 
The (nonparametric) diffuse prior equation A.28 yields; 
1 : 
P -P i=l,...,m(k) 
(A.30) 


y[m)]=! m) °° 
(1-P,-P,) i=0 


293 


The nonparametric model in eguation (A.30) can be obtained from eguation 


A.29 by setting; 


1-20 


70 (4.31) 


i.e., replacing the Poisson parameter with the sample spatial density of the 
validated measurements. The volume V(k) of the elliptical (i.e., Gaussian-based) 


validation region is given in equation (A.5). 


The Parametric PDA: 


Using equations (A.29) and (A.25) with the explicit expression of the 
Gaussian PDF in equation (A.23) yields, after some cancellations, the final 


equations of the parametric PDA with the Poisson clutter model; 


n i=1,...,m(k) 
b+) e ; 
j=1 
Blk) = ; (A.32) 
== i=0 
b+) e ; 
j=l 
where; 
ie gre SGT nA (4.33) 
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1/2 = -P; 


b=1-|27-S(k)| (A.34) 
Pp 
The last expression above can be rewritten as; 
b- ZE) avae A (4.35) 
A Py 


P,: Probability of detection of a target originated measurement, 


P,: Probability of measurements in the gate. 
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